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Abstract 


This report attempts to document the broad scope of issues that must be satisfactorily resolved 
before one can expect to methodically obtain, with a reasonable confidence, a near-optimal robust 
closed loop performance in physical applications. These include elements of signal processing, 
noise identification, system identification, model validation, and uncertainty modeling. Based 
on a recently developed methodology involving a parameterization of all model validating uncer- 
tainty sets for a given LFT structure and noise allowance, a new software, Uncertainty Bound 
Identification (UBID) toolbox, which conveniently executes model validation tests and determine 
uncertainty bounds from data, has been designed and is currently available. This toolbox also 
serves to benchmark the current state-of-the-art in uncertainty bound determination and in turn 
facilitate benchmarking of robust control technology. It is hoped that a convenient new toolbox 
will encourage extensive application of this baseline methodology and toolbox so that researchers 
can build on it to make robust control a truly useful tool. To help clarify the methodology and 
use of the new software, two tutorial examples are provided. The first involves the uncertainty 
characterization of a flexible structure dynamics, and the second example involves a closed loop 
performance validation of a ducted fan based on an uncertainty bound from data. These examples 
also help describe the many factors and assumptions that determine the degree of success in ap- 
plying robust control theory to practical problems. The results demonstrate the highly non-unique 
nature of model validating set of uncertainties, which require only matching predicted to measured 
input and output data, and hence almost any internal model structure is possible. The results 
on the effects of erroneous or non-existent uncertainty structure indicates that model validation 
conditions satisfy easily but often gives unpredictable uncertainty norm bounds, due to differing 
uncertainty structures. 


1 Motivation 

1.1 Robust control theory, a promise 

In robust control theory of the /i variety (for example, see [1], [2]), uncertainties in a math- 
ematical model representation of a physical system are addressed by considering a bounded 
set of models. A widely accepted representation of a bounded set of input-output models 
for general, multivariable, linear, time-invariant (LTI) system is the linear fractional trans- 
formation (LFT) representation 

Q = Fu{Po(W), A u ) Q , Aw := 1FA U G V w 

Figure 1 shows the standard LFT connections of an uncertainty model, A u , augmented 
nominal plant, P Q , and controller, K. The component uncertainty sizes are denoted by the 
block-diagonal uncertainty weight matrix, W, which corresponds in structure to the unity 
norm-bounded block-diagonal structured uncertainty set 

Pb '■= {A u := {diag[5]/ ni , . . . , S r I nr , A 1; . . . , A^] : S t , Aj G IZHoo}, || A^H^ < 1} 

With the above notions of LFT’s and uncertainty structure, a set of plant models can then 
be represented as 

V A := {F U {P 0 (W), A u ), A u G V B } 

and the notion of robustness can be defined such that if a characteristic holds for every plant 
P G V A , then controller K is said to be “robust”. As an example, “robust performance” 
characteristics include both internal stability and a minimum performance level. The amaz- 
ing thing about // robustness theory is that nominal performance (NP), robust stability (RS), 
and robust performance (RP) requirements for the general LFT system shown in Figure 1 


1 





G 


G := F,(?(W) , K) 


Figure 1: Standard LFT uncertainty- plant-controller representation. 


can be succinctly quantified by the following /i conditions: 

NP <=>• NS + || <^22 ||oo <1 ^ d'(G ? 22) = Ma p (G ! 22) < 1, Vcn 

RS NS + F U (G, A) stable, VA e R A 4* n Au (Gu(W)) < 1, Vu 

RP RS + \\F U (G, AjUoo < 1, VA G B A <=> [m a (G(W)) < 1, Vu 


where 


and 


A : = 


hA (G) — 


1 

miriAen {a (A) | det (I — G A) = 0} 



G T>, A p unstructured, complex, performance block 


The point here is that robust control theory is very elegent and computationally reliable [3 
but depends heavy on the uncertainty structure and weight, W. Notice that the assumec 
uncertainty structure and weight is implicit in the augmented nominal plant and its LFT 
interconnectons. 


1.2 Robustness dependence on uncertainty model 

Even after almost two decades of well established theory, there is a persistent gap across 
robust control theory and technology. This means that its usefulness to practical systems is 
limited and is painfully evident from simulation studies of simple toy models to applications 
on more realistic systems, some of which is illustrated in this study. Figure 2 hints at why 
this may be so - a priori knowledge of a physical application is imperfect and there are 
limitations to measurements. 1 

Interestingly the exploitation of a given uncertainty model, with explicit assumptions on 
the system structure, is the key advantage in robust control as well as its main predicament 
in applications. With a given uncertainty model, robustness analysis can address what-if 
questions such as worst case response which can be used as a means to tradeoff system 
robustness with performance. Furthermore, with a given uncertainty model, robust control 
synthesis can build-in robustness with respect to this uncertainty model. In the end however, 
the physical significance of these analysis and synthesis results depends on the relevance of 
the given uncertianty model with respect to the physical system in question. 

1 A plausible historical explanation for this state appears to be that while feedback control was originally motivated 
by physical applications, control theory apparently evolved in its own sphere very often without the distractions caused 
by physical realities and a myriad of peculiar characteristics in physical systems. Modeling is not the primary concern 
of control systems theorists, which makes new control theories more and more difficult to apply. 
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Figure 2: Factors contributing to robust control law performance in real problems. 


These reasons motivate the goal and the particular approach reported here. Namely, we 
seek to improve the usefulness of robust control by constructing uncertainty model bounds 
based on both a priori information and test data, the latter which is usually not used clue 
to the lack of a sensible methodology. 
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2 Structured uncertainty bound determination from data 

2.1 Assumptions on signals and systems 

In dealing with signals for model validation, sufficient care must be given to prevent inad- 
vertent errors. For example, we may not necessarily want to apply tapered time windows for 
input and output signals to improve their spectral spreading and leakage as commonly done 
in spectral analysis, at the expense of spectral distortions in the predicted output. In gen- 
eral, we have to avoid nonperiodicity in the signals because it can lead to significant errors 
when used in predicting output DFT’s, as described in the following subsections. Hence, 
implementing experiments that produce periodicity in the input and output signals appear 
necessary in a proper application of model validation tests based on matching predicted 
output spectra to measurements based on time-limited samples. However, if a periodic time 
sample cannot be attained, the nonperidicity induced errors can be mitigated to some degree 
by applying tapered time windows with a nominal level of spectral distortions. In addition, 
to develop theory with some level of rigor, simplifying assumptions on the system are also 
required, although not more limiting than /i theory in robust control. 


2.1.1 Model of physical system and signals 


Consider a given physical system modeled as a linear transfer function matrix, G, with inputs 
and outputs denoted by u and y respectively and driven by unknown exogenous disturbances 
denoted by d as shown in Figure 3. Our interest is of course the model validation (or 
invalidation) of the system model, G, using a judicious set of inputs and measured outputs 
from an unknown “true” plant model, G trU e ■ Assume discrete time signals so that hereon we 


unknown exogenous disturbance 




reference time with 
unknown initial conditions 


u[m] 




d[m] 


time(m) 


start of test input 


,Ai 2 

end of test input 


Figure 3: Schematic of signals collected for model validation 
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denote, for example, the input u[m] := u(t = mr ) to mean the input at time t = mr where 
r denotes the sampling time and m is an integer index denoting time. While the unknown 
exogenous disturbance signals, d, are everlasting , the test input, u , starts at time k Q and 
ends at k\ . The interval from ki to ko marks the transition from the end of forced response 
to beginning of zero input steady-state response. For more details on signals and systems, 
see for example [4]. 

Suppose the system is modeled as an n-th order state-space, linear, discrete time-invariant, 
and asymptotically stable so that its j'-th eigenvalue denoted by 7 j satisfies 

Assumption 1 : | Tj I < 1 , j — 1, ■ • • , n 
with its output at time m := k Q + k, given by 

n 

y[m ] = E Cj'yJ 1 +u[m ] * g[m\ + d[m ] * h[m\ (1) 

3 = 1 

zero input response 


where {cj,j = 1, . . . ,n} denote constants dependent on generally unknown initial conditions 
at time m = 0, denotes convolution operation, and g[m\, h[m] denote pulse responses 
of models for the physical system and noise filter respectively. Assuming these models are 
causal, 

Assumption 2: g[m] = 0, h[m\ = 0, Vm < 0 

and the convolution terms in equation 1 reduce to, for example 

00 m 

u[m\ * g[m\ := u[l}g[m - Z] = u[l\g[m - l } (2) 

l = — 00 l= — oo 

To deal with the unknown initial conditions, we assume that the corresponding reference 
time for the asymptotically stable system is sufficiently far back relative to time k Q such that 
its zero input response component at time m > k a is insignificant. Specifically, given an e > 0, 
we assume that there exists a positive integer k such that if k Q > k, then j E”=i Cjlf I < d 
Vm > k a . Simply put, 


n 

Assumption 3: c jlT ~ 0) m>k 0 

3 = 1 


For convenience, the reference time point is shifted forward by k Q and hereon for simplicity, 
with assumptions 1 to 3, the output equation 1 is reduced to 

y[k] = u[k] * g[k] + d[k] * h[k), k > 0 (3) 


2.1.2 Model validation using time-limited data 

For model validation in the frequency domain, based on input and output data sequence 
which are of W-points duration (corresponding to a rectangular time window in Figure 3) 

{y} N = {y[ 0], ■ ■ ■ , y[N - 1]}, {u} N = {«[o], . . . , u[n - 1]} 
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consider the ^-transform of equation 3 


Y[z] = G[z]U[z] +H[z]D[z] 


( 4 ) 


where 

OO 

Y[z] := Y. #]*'* 

k=— oo 

and similarly for U[z) and D[z] while G[z\ and H[z] denote transfer function matrices whose 
frequency response functions are defined on z = e lWT . With time-limited data, we cannot 
compute the necessary ^-transforms in equation 4 which motivates us to rewrite as 

— 1 TV— 1 oo 

y[z) = Y y\ k \z~ k + E yl k \z~ k + E y^ z ~ k 

k=— oo k= 0 k=N 

t — 1 TV— 1 oo \ 

u[k]z~ k + u[k]z~ k + u[k]z~ k I + H[z]D[z] 

k=— oo k= 0 k=N / 


The above relation rewritten only at discrete set of frequency points evenly spaced on the 
unit circle in the £-plane is 

| . 27 r 

Y N [z r \ — G[z r ]U N [z r \ + RN[z r \ + — 7 =H[z r \D[z r \, z r = e l ^ r , r = 0, . . . , TV - 1 (5) 

v iV 

where Discrete Fourier Transforms of samples { y} N and are given by 2 


N - 1 


iV-1 






E 


,-fc 


Ui\r\z r 1 : — 




E “W ; 


,-fe 


and the mostly unknown residual 

1 f ^ 


-^/vl^rl • 


v/N [ 


<1 - Y y ^ z r k - J2 y ^ z r k + u ^ z r k + U 

k=N 


,-k 


k =— oo 


\k =— oo 


fc=7V 


is due to time-limited input-output data. Ignoring for the moment the effect of unknown 
exogenous disturbances and noise, -^i7[z r ]D[z r ], an upper bound on the residual can be 
written in terms of the worst case input signal and a measure of system gain as [5] 

2 00 
\R N [z r \\ < -7=|M|oo 'Yj 

k=Q 

Since the DFT of an TV-point sample is equivalent to the Fourier Transform of the corre- 
sponding everlasting TV-periodic signal, this means that the residue consists of the difference 
between an TV-periodic output constructed from measured output, and the calculated output 
from a given discrete-time filter driven by an TV-periodic input constructed from measured 
input. Of course without a priori knowledge of the given measured signals, there is no reason 

2 In MATLAB-Signal Processing Toolbox, the scale factor i is used in the DFT calculation (rather than 
implied by our DFT and inverse DFT definitions) and a unity scale factor in the inverse DFT expression. 
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to assume that the constructed iV-periodic sample is representative of a nonexistent peri- 
odic input-output response of the physical system. Also, this residual holds for a generally 
unavailable “true” model (the given discrete-time filter) which produced the sample time 
records, hence, additional model error terms must be included subsequently for a model 
validation test. 

If the measured input and output signals are truly iV-periodic, then 

Assumption 4: RnI^v] = 0, Vr 


A time-limited signal, for instance {y} N , can be seen as the product of everlasting time 
signals, y, with a rectangular window function, w of width N, i.e. 

VN[k\ ■= y[k]w[k], k = — oo, . . . , oo; Y N [z] ■= Y[z] * W[z] 

where {t/jv} := {• ■ ■ , 0, {y}^, 0, . . . , } is the everlasting time signal which contains the time- 
limited signal. In spectral estimation problems tapered time windows are often used to 
mitigate spectral spreading and leakage caused by the rectangular window function which 
can be seen as a convolution with a sine spectrum. Basically, the tapering works because 
the signal spectra in some bandwidth of interest is minimally distorted unlike at other 
frequencies. 

In contrast to the spectral estimation/analysis problem, we are interested in accurately 
characterizing the predicted output error defined by (see Figure 3) 


E[z r ] - Y N [z r } - Y n [z t ], r = 0, . . . , JV - 1 


( 6 ) 


where 



Yjsf \z r \ 


Z{y[k]w[k\) = Y[z r ] * W[z r ] = (G[z r ]U[z r \ + Noise) * W[z r ] 

G[z r \Z(u[k]w[k]) H — — H[z r ]D N [z r \ = G[z r \(U[z r ] * VF[z r ]) H — — H[z r \D N [z r 
VN v N 


The term Yn[z t ] denotes the ^-transform, denoted by Z(-), of windowed measured outputs 
evaluated at z r while the term Y N [zk\ denotes the predicted outputs based on assumption 
4, i.e. R.N[z r ] = 0 in equation 5. The predicted output error given by equation (6) with 
TV-periodic assumption reduces to 


E[z r ] = (G[z r }U[z r ]) * W[z r \ — G[z r \{U[z r \ * W[^ r ]) + Noise Terms, 


r — 0, . . . , N — 1 (7) 


so that even for the true model, i.e., G[z\ = G[z\, and excluding noise effects, the above 
predicted output error will not be zero; only a rectangular window will indicate a zero 
predicted output error for the true model. Intuitively, it matters whether the time windowing 
is applied at the input to the plant or at its output. This means that tapered windows 
should be used sparingly in model validation for preconditioning signals (for example to 
force periodicity or even mitigate aliasing). 

In the context of conducting model validation experiments, we can try to satisfy the 
A^-periodic assumption by either (i) literally implementing an .A-periodic test input or (ii) 
implementing a test input sequence that will produce a rest-to-rest Appoint sample, for 
instance by beginning an Appoint sample such that the zero input response from earlier time 
is negligible and by including at the tail end of the sample the decaying output response 
(due to zero input) until steady state. In either case, the experiment must be designed such 
as to allow a certain prescribed periodicity in a measurement sample, without using tapered 
windows. The effects of unknown exogenous signal terms are difficult to account for and is 
dealt with by prescribing noise/disturbance allowances in the periodogram of the appropriate 
signals based on a priori knowledge of say their estimated spectra. 
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2.1.3 Example 


In this section we illustrate possible errors in the predicted output even if the model to be 
validated is the true model which generated the signals. Suppose the true system is a lightly- 
damped second-order discrete time plant defined by G[z] = zp2sys( [— 1, — 1], [.9 + j *.4, .9 — 
j * .4], 1) with sampling time r = .1 sec. Figure 4 shows the Bode plot of the true plant and 
simulated input and output time records imbedded with a low level of additive output noise. 
The test input chosen is simply a random number sequence having a uniform distribution, 


G[z] = zp2sys([-1,-1],[.9+j*.4,.9-j*.4],1) 



Frequency (radians/sec) 


Sampling Time =0.1 

4 1 r 1 1 1 1 r 



0 50 100 150 200 250 300 350 400 450 500 



output noise (-.), o = 0.0001 
150 1 1 1 1 1 1 1 



0 50 100 150 200 250 300 350 400 450 500 

Time (sec) 


Figure 4: Frequency response of G[z] (left), simulated input/output data with low noise (right). 

zero-mean with unit standard deviation generated by inp = [zeros(1000, 1); randn(l000, 1); 
zeros(l000, 1); randn(l000, 1); zeros(1000, 1)], independently for each case. Similarly, the 
assumed output noise considered are generated by uniformly random time sequences having 
zero mean and with two different levels of standard deviations, cr = 10 -4 to represent low 
noise case and a = 1CT 1 for high noise cases. Starting at rest, it is assumed that the input 
is periodic over exactly 2000 points (or 200 seconds) of which the latter half consists of zero 
input to let the system decay to near zero state before the next period begins. Of course 
in principle, it does not matter where the 2000 points begins and ends for periodicity but it 
must be exactly 2000 points, else, the errors can be huge. For this reason, in applications, 
it helps to select the end timepoints of a sample such that the measured output do not end 
abruptly at both ends. 

Table 1 show cases under which the predicted output errors are computed with variations 
in output noise levels, signal periodicity, and window function. In all cases, the true model, 
which was used in generating the simulated signals, is used in computing the predicted output 
errors with only additive output noise - obviously to examine potential hazards without even 
worrying about model or parametric errors. 

The top two figures in Figure 5 shows the input and output time records and the cor- 
responding predicted output error for low noise, periodic Cases LPR and LPH. The only 
difference between these two cases is that the latter case uses a tapered time window. In 
both cases the time samples are chosen to be perfectly periodic over 2000 points, starting 
from point 1001 and ending with point 3000. A sufficient time interval of 100 seconds with 
zero input occurs at the latter half of these samples to satisfy periodicity with negligible 
zero input response from the past. The third subfigure for each case show the spectra mag- 
nitudes | F/v [zr] | (solid line), \G z r ]I7jv[zr]| (dash-dot line), |Y/v[z r ] — G[z r ]U n[z t \\ (dash line), 
and | -Dj\r[z r ] | (dotted line) while the noise filter H[z r ] is assumed to be unity. In all cases 


Case 

Noise Level, a 

Periodicity (point range) 

Time Window 

LPR 

1CT 4 

Periodic(1001-3000) 

Rectangular 

LPH 

10" 4 

Periodic(1001-3000) 

Hanning 

LNR 

1CT 4 

Non-periodic(1001-3700) 

Rectangular 

LNH 

1CT 4 

Non-periodic(1001-3700) 

Hanning 

HPR 

HT 1 

Periodic(1001-3000) 

Rectangular 

HPH 

1CT 1 

Periodic(1001-3000) 

Hanning 

HNR 

1CT 1 

Non-periodic(1001-3700) 

Rectangular 

HNH 

1CT 1 

Non-periodic(1001-3700) 

Hanning 


Table 1: Cases simulated for predicted output errors. 


considered, the spectra magnitudes Rv[^r] and |G[z r ]Hjv[£ r ]| appears to overlap in the log- 
arithmic scaled figures and yet their corresponding predicted output errors (dash line) are 
not insignificant. 

The spectra magnitude plots for Case LPR shows that the predicted output error is to- 
tally due the output noise (the dash and the dotted lines overlap). However, if the periodic 
samples are Hanning windowed as in Case LPH, the predicted output error increases signif- 
icantly relative to the output noise spectra floor. Furthermore, as indicated in the previous 
figures the frequency shape of this increase is dependent on the plant transfer function, i.e. 
dependent on the frequency convolution of Hanning function with plant transfer function. 

The bottom two subfigures in Figure 5 corresponds to Cases LNR and LNH which illus- 
trates the consequence of violating periodicity assumptions in the signals. This particular 
nonperiodic time samples start from point 1001 and ends at point 3700. The predicted 
output error for Case LNR is due to the nonperiodicity of the signals used in computing 
the DFT’s and is significantly larger than the predicted output error in Case LPH which is 
primarily due to tapered window at the input for a periodic signal. With the application of 
the Hanning window to the nonperiodic signals in Case LNH, the predicted output error is 
significantly reduced when compared to the rectangular windowed nonperiodic Case LNR, 
but still the error peak near 4 rad/sec did not improve. Interestingly, the predicted output 
error level in Case LNH is comparable to Case LPH, whose errors are due to the generally 
nonequivalent frequency convolution at the input and the output for a general window func- 
tion. This is likely due to the fact that applying the Hanning window on the nonperiodic 
signals in Case LNH makes it somewhat periodic, but still the error due to nonequivalence 
of the frequency convolution at the input and output remains, similar to Case LPH. 

To illustrate the effect of high output noise level, Figure 6 show Cases HPR, HPH, HNR, 
and HNH. Similar to the ideal Case LPR but with a higher output noise floor, Case HPR 
shows that the predicted output error equals the increased noise floor. The predicted output 
errors in Case HPH shows the combined effects of the input/output Hanning windowing 
(comparable to Case LPH) and the higher noise floor (comparable to Case HPR) so that 
at low frequencies, the noise floor dominates the predicted output error while near plant 
resonance the Hanning window induced frequency convolution error dominates. The effects 
of using nonperiodic signals in Cases HNR and HNH is analogous to the low noise cases, 
namely, a large predicted output error is incurred by the nonperiodicity (Case HNR) and 
is subsequently improved by Hanning windowing (Case HNH) - of course subject to the 
limitation imposed by the higher output noise floor. 

In summary, this example demonstrates that even with no model errors, significant errors 
in the predicted output can result due to the effects of nonperiodicity, windowing, and noise. 
If the signals are known to be nearly periodic, rectangular windowing appears to reduce the 
predicted output errors, however with nonperiodic signals including any signals with high 
known levels of noise, as in the Ducted fan example to be discussed later, tapered windowing 
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Figure 5: Cases LPR (top left), LPH (top right), LNR (bottom left), LNH (bottom right): In- 
put/output records and predicted output error. 


appears more suitable. Of course in actual applications, the added complication of model 
errors is inevitable. So, in model validation for robust control, we utilize the uncertainty 
structure and plant set associated with robust control to match the measured output with 
any one predicted output amongst a given set of models. 
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2.2 Model validation framework 


We begin with making Assumptions 1 to 4 to restrict the class of signals and systems under 
consideration. Motivated by multivariable robust control theory (see for example [1]), we also 
limit the model structure to those that can be represented by LFT connections. The basic 
model validation question involves determining whether a particular plant model, belonging 
to a given set of plant models, and a particular exogenous signal, belonging to an unknown 
but bounded set of exogenous signals, could have produced a given pair of measured input and 
output signals. Going beyond this basic question, we also show how reasonable uncertainty 
bounds can be constructed for use in robust control or performance validation. 

Figure 7 shows the basic LFT framework and its canonical form for model validation. 
Notice that for a closed loop system model validation where r and K are assumed known, 


unknown exogenous disturbance 



i>=m 

r 


Figure 7: Model validation framework and canonical form. 


e u := u — u = 0 if e y = 0 so that only the output error need to be considered. Given mea- 
surements of the output, y, command input, r, an augmented nominal plant and disturbance 
filter, P, controller, K, output noise filter, V, and a set of bounded structured uncertainty, 
T>w, the set of plant (robust control design models) 

Vw '■= {F U {P, A), A G XV} 


is said to be a model validating set if at each frequency (uj k = k = 0, . . . , N—l) it contains 
an uncertainty model A G XV such that there exists exogenous disturbance signals, e and 
v with 1 1/? 1 1 < 1 for which 


y = F u (G(P,K,V), A) 


e 

v 

r 


and 


V ■= (A G C mxn : A = diag(h 1 X mi , . . . , 5 r I mr , A r+1 , . . . , A r ), G F h A k G C mM } 
W := diag Wi, . . . ,w T X T ) 

XV := (A G V : A = W A B ) 

A B ■= {AGP: b(A) < 1} 


With Assumptions 1 to 4, the set of signals y, r, r), u, are viewed as the DFT’s of their re- 
spective time-limited samples while systems, P, K, V, G, W, are viewed as transfer function 
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matrices evaluated at the same uniformly spaced frequencies np to Nyquist frequency. The 
symbol T u denotes an upper LFT structure and G(P, K, V ) denotes the augmented closed 
loop transfer function matrix containing both noise and disturbance filters. Notice that we 
do not require A G VSHoo so that stability and causality of the perturbation are not required, 
only that it is norm bounded. The terms 5j, i = 1, . . . , r, denote repeated scalar uncertainties 
which in general can depend on frequency in contrast to parametric uncertainties which are 
implicitly assumed to be unknown constants. 

This question was first posed in [7] and they have shown that the validation test can be 
formulated analogous to a skewed-p or alternately as a quadratic optimization problem. For 
the “most commonly applied cases: additive and input or output multiplicative” [8] it has 
been shown (and in time domain [9]) that the validation test takes the form of a convex 
feasibility problem when the fictitious signal 7] do not depend on £. More recently, [10] has 
shown that by choosing a plant uncertainty model based on coprime factorizations of the 
nominal plant and known controller, all closed loop transfer functions can be made affine in 
the coprime factor uncertainties so that the convexity of the validation test can be preserved 
through enforcing G(P, K, V)n = 0. In dealing with closed loop validation in [10], the noise 
model enters as an additive output outside the loop and hence only Gn is relevant. On the 
other hand, in the work reported in [9], an open loop model validation problem with additive 
output noise is considered so that only Pu is relevant. However, it is not clear what the level 
of conservatism is on this particular form of plant uncertainty structure especially since no 
additional user specified uncertainty structure can be accomodated nor is necessary. 

A key problem arises if one wishes to test a more general uncertainty structure for model 
validation beyond the “most commonly applied cases: additive and input or output multi- 
plicative” (as mentioned in [8]), namely, the loss of convexity in the feasibility problem. This 
means that from a numerical implementation standpoint, the validation test is effectively a 
sufficient condition. Furthermore, if the validation test is passed for a particular plant set, 
there exists an infinity of other sets which will also be validating, irrespective of its convex- 
ity. In our view, this non-uniqueness diminishes the significance of a validation result for a 
particular set. 

For these reasons, we take an alternate view of the model validation problem, as described 
in more detail in the next section. Instead of trying to validate a specific set with respect 
to given data by asking “Is T>w model (in) validating?” we ask instead “Does a model 
validating set T>w exist?” The former question is a test on a specific uncertainty set resulting 
in a positive or negative answer while the latter question is a test to determine whether 
some uncertainty set with the given LFT structure exists that can satisfy model validation 
conditions. The latter question is of course a necessary condition to the former and it turns 
out to be much simpler to test. Furthermore, if a finite size model validating set T>w exist, 
then it turns out that all model validating sets having the same LFT uncertainty structure 
can be readily parameterized [6]. 


2.3 3-Step Approach to model validation and parameterization 

A new approach breaks up the model validation question into 3 parts: 3 

1. Does a pair of signals (£,/3) where ||/3|| < 1 exists such that e y = 0 ? 

2. Parameterize the set of signals 

'■= {[£(^),/3(0,VO>^(0>VO] : e y = 0, </> G $,-0 G T} 

3 For more details see [6] 
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3. Does a A G D exists such that £ = A?/ for some signal (£,?/(£, (3)) G ? 


Figure 8 illustrates the question posed in Step 1. Such a pair of signals exists if and only if 



Figure 8: Does (£,/?) exists, where ||/3|| < 1, such that e y = 0 ? 


e° y G Im(Af) \ . . 

||7f (M+)^ e °|| < 1 J Va;fc (8) 

where the constant matrices are given by 

&y • Vmeas G 2 3T" 

M := [G 2 i,G 22 ] 

Im(N M ) = Ker(M) 

T 2 := Im{{N M ) v ) L 

Notice that M is a crucial matrix which captures the uncertainty freedoms to be used for 
model validation. This test involves only a constant matrix check at discrete frequency 
points of interest in contrast to a (non)convex feasibility problem. If the above necessary 
condition cannot be satisfied, this means that no matter how large the uncertainty radii W 
(imbedded in G) is, the model is invalidated because either the a priori uncertainty structure 
is too restrictive or the noise allowance is too small. Hence to proceed in trying to construct 
a model validating set, one has to improve the nominal plant model, P, and/or modify 
the uncertainty structure T> to increase its richness, and/or increase the noise allowance by 
modifying the noise filter model V. 

If the above constant matrix test passes, then one can proceed to Step 2 and parameterize 
the set of signals 

■= {[^(0,'0),^(0,'0),/5(0,'0)] : e y = 0, 0 G G T} 
where all feasible triples (£, e, v) are given by 



where 0 is arbitrary and 0 satisfies 


and 


IMI < K := 0 - l|T"(M+)„eJ|| 2 


( 9 ) 


ft 


0 
u 

( Nm)/3 


N m U 


S/ 1 0 

0 4,, 


[M + — N m ((N M ) p) + (M + ) p\ e° y 


[4 t 2 ] 


Si 0 

0 0 


U H 


( 10 ) 

( 11 ) 

( 12 ) 
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Figure 9 illustrates the final step of the existence test, which is to check whether a A E V 
exists such that £ = Arj for some signal pair (£, ■/?(£, (3)) G S ( p^ where rj = Gn£ + Guf3 + 
Gi^r. ft turns out that this part of the existence test is satisfied if and only if there exists 


tK<IW) 



^(qw) 


Figure 9: Final step of existence test 

(£,0) G Sfaf such that £i = 0 or r/j(£,/5) ^ 0, Vi, i.e. ©-realizable, and for each repeated 
uncertainty block, (£*, Vi{£i 0)) are collinear. For the special class of uncertainty structure 
with no repeated uncertainties, i.e. r = 0 but not requiring G(P, K,V)n = 0, with a 
satisfaction of the necessary conditions for existence in equation 8, it is almost certain that 
there exists (£,/?) G such that rji(£,/3) ^ 0 , Vi. This means that for this important 
special class of uncertainty structure, the necessary conditions in equation 8 are actually 
necessary and almost sufficient for a model validating set V w to exist. 

If a model validating set Vw exists, then all model validating sets of plants are given by 

'Pwcfnp ■= {©«(©, A), A G V w } 

where 0 G C 0 G C n 0 ||0|| < b a , and W := diag(wi/ ni , . . . , w T I nr ) is any matrix satisfying 



116(0,0)11/, | 
11^(0,0)11 th 

, i = 1, 

...,r 


(£, v) is ©-realizable and 

are parameterized by 





& £>o,i ^ i ^ 

* \ 
0 / 


(13) 


Vi = Vo,i + [G 11 

G^n | 


(14) 

dist ( “ Fi) (£ i , rji) = 0,i = 1, . 

. . , r, and 





Vo [f?n G 12 ] | j” + Gi^r. (15) 

It is intuitively pleasing to note that if model validating sets exists, they are highly non- 
unique. 


2.4 Validation as a feasibility problem 


Using the previous parameterization, we re-examine the earlier question, “Is V w model 
(in) validating?” Suppose a model validating set V w exists. Since all model validating sets 
for the chosen uncertainty structure are given by the above parameterization, the original 
question can be posed as follows: 

Given W, at each frequency does a (0,0), ||0|| < b a exists such that 


116 ( 0 ; 0)11 < | 

11 ^( 0 , 0)11 - ' ’ 


i = 1, ■ • • ,r 


where (£, r/) is ©-realizable and dist^^(£j, Vi) — 0, i — 1, . . . , r ? 


(16) 
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The T >- realizable and collincarity conditions can be implicitly satisfied by the following form: 
Given W, at each frequency does a (0, 0), ||0|| < b Q exists such that 
6(0,0) <- (0,0) = 0, 5 t <E T % 

| 6 | - \wi\ < 0 

116(0, 0)|| 2 - M 2 ||?7i(0,'0)|| 2 < o, i — r + 1, . . . , r ? (18) 

Although 6 and r\i are affine in 0 and 0 as given in equations 13 and 14, the feasibility con- 
ditions are in general nonconvex. Specifically, the collinearity requirement due to repeated 
(and/or real) uncertainty involves a quadratic equality since 6 is unknown equation 17 and 
the norm inequalities for the nonrepeated uncertainties although quadratic involves a differ- 
ence in norms (equation 18). Nevertheless, one is hopeful towards a workable methodology 
if one observes from equations 17 and 18 that a feasible solution can be generated by first 
focusing on finding a pair (0,0) which satisfies the collinearity condition in equation 17 and 
then easily selecting sufficiently large weights wy, i = 1 , . . . ,r, to satisfy equations 17 and 
18. 



2.5 Optimizing for smallest scaled set 

A useful set of model validating plants must be “small but not smaller”, i.e. , it should not 
be an unnecessarily large set of plants which will limit performance but it must contain 
the “true” plant. Given a candidate set XV of uncertainty norm radii, W, we seek a 
corresponding smallest scaled set T> x w which is model validating. This question can be 
posed by imbedding the previous feasibility problem in the following optimization: 




min 

X 2 




i,...,5 r ,x 2 



subject to 

6(0,0) 

- 6^(0,0) = o, 6 

6 2 — x 2 \wi 2 < 0 

: e Ti \ . 1 

> « = 1 ,...,r 

(19) 


116(0 , 0)ir 

-* 2 H 2 lk(0,0)ll 2 

< 0, i — r + l,...,r 

(20) 



II6II 2 

< b 2 a 

(21) 


Although the above problem is not convex in general, the conditions involve polynomi- 
nals in 0, 0, Si, . . . , S r , x' 2 with are at most cubic order. The existence of a feasible point 
(0, 0, Si, . . . , S r , x 2 ) in the above optimization problem is equivalent to the existence of a 
model validating set V x w ■ Of course when x < 1, the original candidate uncertainty set, 
XV satisfies model validation conditions. Notice that if there are no repeated uncertainty 
(r = 0), any choice of 0 and 0 where ||0|| < b Q will likely admit a model validating set 
T> x w, since a larger uncertainty radii, W will admit a feasible point. However, even without 
repeated uncertainties, the constraints are in general not convex in 0 and 0 and therefore 
the numerical optimization task is likely nontrivial. 

The important issue of how one would simultaneously choose apriori weights for paramet- 
ric (uii, i = 1 , ... ,r) and nonparametric (wj, i — r + 1 , . . . , r) uncertainties to form candidate 
uncertainty weight is beyond the scope of this work. Whatever is chosen in an application, 
these candidate uncertainty weights should reflect the relative importance of the uncertainty 
components by design, based on additional a priori information on the particular application. 
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2.6 Convex model validation problems 


For the special case when [Gn, Gi 2 ] = 0, 77, becomes independent of £ and 0 (or equivalently 
0 and 0) so that the term ^(0,0) simplifies to a constant 77* = Gi 3i r for given data and a 
convex feasibility/optimization problem results. This means that in general, both the choice 
of the LFT uncertainty structure and the way unknown exogenous signals affect the model 
determines whether a convex problem arises. For this special case, the problem of finding a 
smallest model validating unmodeled dynamics reduces to: 


subject to 


min x 2 

<t>,4>,5 i,...,S r ,x 2 

116(0, 0)|| 2 -x 2 |^| 2 ||G 13i r|| 2 < 0, 

IHI 2 < bl 


i = l,...,r 


and since 6(0, 0) is affine and the inequalities represent ordinary norm bounds, a feasible 
set if it exists will be convex. For computation, we rewrite these as an optimization problem 
involving a linear cost function subject to a set of linear matrix inequality constraints [11]: 


min c T z 

Z 


subject to 


Qi? - 116, ill 2 

sym 

S iZ 

I 

r bi 

sym 


Lz I 


> 0, i — 1 , . . . , r 

> 0 


where 2 : = 


[Re(0); lm(0); Re(0); lm(0); x 2 ] E fl 2n ^ +2n <t>+ 1 and 


c 

Qi 

Si 

L 

A 


[0(2^ +2 n f y, 1] 
r 2Re(6"0 ?: kl), H 2 ' 
Re(h2j) — Im(Oj] 
Im(Oj) Re(h2j) 


IGiaTlI 2 ] 


T Re (A) 

,0 


Im(A) 



0 


2nfk x 2 n. 


"ijj 


‘2n„s 


0 2 


2n^x 1 


0 0 jln^ 

I fty J Si,:, 0 0 


( 22 ) 


In a recent study, [10] notes that even if a selected uncertainty structure in open loop is made 
to satisfy [Pn, P12] = 0 (see Figure 7), the presence of feedback can lead to [Gn,Gi 2 ] 7^ 0 
resulting in a nonconvex feasibility test. In general, it can be observed that since 7] = 
[Gn, Gi 2 ](|) + Gi 3 r, [Gn, G i2 ] = 0, if and only if all of the following conditions hold: 


1. each uncertainty block is not connected to the plant in a feedback, 

2. uncertainty blocks are not connected in such a way that an output signal from one 
block could influence another block, 

3. all exogenous input signals do not influence any uncertainty block 


These conditions are for open loop plants so that for systems with embedded feedback 
controllers forming an augmented plant, some of the above conditions may not hold. Clearly 
these conditions, which are sufficient conditions for convexity, can be very limiting in the 
choice of the uncertainty structure and exogenous input connectivity. A convenient to use 
interior-point algorithm for convex problems from LMI Control Toolbox [20] can be used to 
solve for the global minimum when a feasible set exists. 
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2.7 Smallest unmodeled dynamics subject to parametric uncertainties 

In the previous section, we dealt with repeated scalar and/or real uncertainties that may 
depend on frequecies. In this section, we consider “parametric uncertainties” which by 
physical reasons are considered to be constants but unknown. One approach to handling 
this unknown constants is to consider it as an independent parameter identification problem 
in the context of system identification, a la [12, 13, 14, 5, 15]. This of course means avoiding 
model validation and LFT uncertainties entirely. The alternative viewpoint we take uses 
these unknown constants as additional variables in minimizing the norms of model validating 
LFT uncertainties. 

If parameters imbedded in the plant changes (as in parameter identification problem) such 
that the crucial model validation design freedom matrix, M = [G 2 i,G 22 ] changes, this can 
cause cascading changes which are difficult to evaluate in a numerical optimization context. 
Hence, we try to avoid changing the augmented plant directly by viewing the imbedded plant 
changes or parametric uncertainties as LFT uncertianties which are allowed to vary without 
changing the augmented plant. 

In the case where we have competing unmodeled dynamics and parametric uncertainties, 
we suggest fixing an allowance for parametric uncertainties while minimizing non-parametric 
uncertainties (of course in addition to a fixed noise allowance). This may be more reasonable 
physically than the approach whereby a priori candidate uncertainty weights are fixed for 
both parametric and non-parametric uncertianties and then determining a smallest scaled 
set. The former approach amounts to a slight modification of the optimization problem for a 
smallest scaled set in the previous section and choosing the x scale in equation (19) to unity. 


2.8 Unknown but bounded exogenous signals 

To deal with unknown exogenous signals such as process and measurement noise, we consider 
them as unknown but bounded. No assumptions on the statistics and the independence of 
the unknown signals are required. In trying to satisfy the model validation conditions, 
these signals are treated as allowances to be used in finding the minimal norm LFT model 
uncertainty necessary. As such, erronously assuming an overly conservative level of noise will 
likely lead to optimistic levels (smaller than actual) of model validating LFT uncertainty. 
On the other hand, assuming a lower level of noise allowance (than actual) will likely lead to 
pessimistic levels of (larger than true value) model validating LFT uncertainty. Consequently, 
it is important to specify a reasonably accurate model of the noise and disturbances in the 
system for model validation and uncertainty bound determination purposes. 

For convenience, the combined noise/disturbance allowance in the canonical form is given 
by 

II/?MI|2<1Vu; ^ m\ lx := II V/3IL < i (23) 

denotes the filter matrix for disturbance and noise allowances and 

0 denotes the exogenous signal as it enters the loop. Specifically, Vp is designed to reflect 

anticipated spectra of 0 and to normalize the unknown norm bound of 0 to unity as in 
equation (23). The assumed block diagonality in Vp implies the independence of e and v. 

In some systems, reasonably accurate models for V e and V v may be available. In other 
systems where typical spectrums of i(u>) and 0(u) are available and if their individual chan- 
nels are known to be independent, a stable discrete time filter can usually be fitted for each 
channel and realized as a diagonal filter. In some open loop unstable systems (such as the 
Caltech’s Ducted fan in a certain flight configuration), a stabilizing feedback controller is 


where Vp := 
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necessary to collect data. Under such conditions, developing suitable models to represent 
unknown noise and disturbance effects from data can be tricky because of the amplifying 
and correlation inducing effects of feedback which is partly determined by loop dynamics, 
which itself is not precisely known. 


2.9 Performance validation using uncertainty model 

A successful performance validation of a control law on a physical system is typically an 
important end goal of any control law design. It is also a definitive test on the validity 
of the mathematical model used in the analysis of the control system, using a single plant 
model or a set of plant models. A key premise in robust control is that any uncertainties 
in the mathematical model is to be contained by an appropriate uncertainty ball around a 
best (nominal) model. With this premise, a robust controller is designed, which is effectively 
an optimal control with respect to a particular uncertainty model. Specifically, a robust 
control law is usually designed by optimizing the predicted worst case performance over 
the uncertainty model. So in practice, it is crucial that the optimized predicted worst case 
performance will actually bound the performance based on measurements. 

The predicted worst case performance over a set of plants defined by an uncertainty model 
over a given nominal is illustrated in Figure 10 where skewed-/! is defined at each frequency 
as [1] 




| minp : det 



0 


0 

piper f 


F t (P, K)AW 


0,cr(A) < 1 


-1 


max a 

unc )<i 


WperfFu^F^P, A'), A unc W unc ) 


(24) 


where A : = 


0 


0 

A perf 


and W 


W unc 0 

0 lUpefj 


In this study, performance 



Figure 10: Predicted worst case performance based on identified uncertainty model 

validation will mean comparing the predicted performance based on an uncertainty model 
to the measured performance. More specifically, the measured performance will be defined 
as the measured maximum amplification over bounded multivariable signals, whereas the 
predicted performance is taken to mean the predicted maximum amplification over a set of 
plants. Since it is a bit difficult to test and measure the response of all possible inputs in 
the hope of directly measuring the worst case signal amplifications, we identify an empirical 
model of the closed loop system from r to e, T er , with the goal of computing the maximum 
singular value of its frequency response matrix, W per fT er . 

In general, any stabilizing control law’s performance could be validated using the same 
uncertainty model, since a useful uncertainty model of a plant should be independent of the 
control law. Indeed, a consistently successful performance validation of several stable control 
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law is a reasonable indication of the validity of an uncertainty model in question. Finally, 
notice that it is a bit difficult to experimentally determine the worst case responses with 
respect to unknown exogenous inputs. Hence, in the above suggested validation, we compare 
the measured and identified worst case responses from command input, r, to measured 
output, e, to the predicted performance n s from r to e only. 


20 



3 Uncertainty Bound Identification (UBID) Toolbox 


3.1 Prerequisites for UBID Toolbox Version 0.1 

• MATLAB 6.0 data structure required 

• Optimization, LMI, and /i Toolboxes required 


3.2 Summary of Commands 

Main routines: 


Command 

congrad 

f ungrad 

mnmvcl 

noise_allow 

ovbndunc 

sseigunc 

ubid 

ubidLLmi 

ubid_uncstr 

weightunc 


Description 

Defines constraint function and its gradient for smallest set optimization 
Defines objective function and its gradient for smallest set optimization 
A smallest model validating set algorithm 
Defines noise and disturbance allowances in model validation 
Interactively fits overbound of computed model validating uncertainty norm 
Forms augmented state space plant with eigenvalue uncertainties 
Sample main program for computing a smallest model validating set 
Sets up convex program of a smallest set optimization problem 
Forms uncertainty structure 

Define normalizing uncertainty weights for optimization 


Utility routines: 


Command 

blkd2x2 
f orm_sample 
pinv_frd 
svd_frd 


Description 

Transforms state matrix into 2x2 block diagonal form 

Divide a time record into windowed-overlapped-zero appended samples 

Computes pseudoinverse of a FRD model 

Computes SVD of a FRD model 


Routines from other toolboxes: 


Command 

fmincon 
f easp 
mincx 
drawmag 


Description 

Solve NLP using an SQP algorithm from Optimization Toolbox 

Solve convex feasibility problem from LMI Toolbox 

Solve convex program from LMI Toolbox 

Interactively fit a specified order stable filter from /z-Toolbox 
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Computational flow of UBID Toolbox 



overbound interactive fit 

uncertainty norm 
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3.3 Command Reference 


3.3.1 blkd2x2 

blkd2x2 


Purpose 

Transforms a state space system into 2x2 block diagonal form. 


Synopsis 

[A,B,C,D,P,nreal] = blkd2x2(a,b,c,d,Ts) 


Description 

The inputs (a, 6, C, d ) are state space matrices and the outputs (A, B , C, D) are block- 
diagonalized state matrices computed from 

A = P 1 * a * P] B = P 1 * b] C = c * P] D = d 

where P denotes a similarity transformation matrix. The output nreal denotes the number of 
pure real eigenvalue/eigenvectors which are grouped in the first nreal blocks followed by 2 by 
2 blocks in the A matrix corresponding to complex conjugate pairs of eigenvalues. These blocks 
are arranged in an ascending order of thier absolute frequencies. For a discrete time system with 
sampling interval Ts as an (optional) fifth input argument, an equivalent s-domain frequency is 
computed from log z. 

-L s 
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3.3.2 congrad 

congrad 


Purpose 

Define constraint functions and their gradients for a smallest set optimization problem. 

Synopsis 

[Cin,Ceq,Gin, Geq] = congrad (x,xio , etao , EF , PEFGH , Wunc ,bo ,blkf , cdimM,rankN , . . . 
rankM , nrrs , nrcs) 


Description 


This subroutine computes the constraint functions and their gradients for a smallest set 
optimization problem which is used by the optimizer fmincon in MATLAB Optimization 
Toolbox. This optimizer uses a sequential quadratic programming algorithm to minimize a 
constrained nonlinear multivariable function. The constraint functions are assumed to be in 
the following forms: Cin(x) < 0 and Ceq(x) = 0. Their corresponding constraint function 
gradients are assumed to be: Gin(x) := d(Cin(x) ) /dx and Geq(x) := d(Ceq(x) ) /dx. The 
output variables Cin, Ceq, Gin and Geq are evaluated at the value X where 


x(l : n ip) 

x(n if, + 1 : 2 n^) 

x{2 + 1 : 2 + n^) 

x(2 + 71$ + 1 : m) 

x(m + 1 : m + n rrs ) 

x(m + n rrs + 1 : m + n rrs 


2n rcs 'j 


x[m + n rrs + 2 n rcs + 1) 


Re(i/j) 

(jl^ X 1) real 

Im( ifr) 

(n^ X 1) real 

Re{(p) 

X 1) real 

Im{4>) 

X 1) real 

dl, ■■■: dn rrs 

[n r rs X 1) real 

Re(dn rra +i), Im(d nrra+1 ), . . . 


R^{,d nrrs -\-n rcs ) 1 d m {d nrrs + nrcs ) 

(2n rcs X 1) real 

x 2 

real 


where 777 := 2 + 271^. The remaining parameters which are defined at a given frequency 
are as follows: 


xio 

etao 

EF 

PEFGH 

Wunc 

bo 

blkf 

cdiM 

rankN 

rankM 

nrrs 

nrcs 


= £o 

= Vo 

= n 

= [Gn, Gi2]S2 
= normalizing unc wts 
= noise uncertainty radius 
= uncertainty structure 
= number of columns of M := [G 21, G22] 

= minimum rank (over freq) of varying matrix (Nm)(3 
= minimum rank (over freq) of varying matrix M 
= number of Repeated Real Scalar uncertainty blocks 
= number of Repeated Complex Scalar uncertainty blocks 


(771 X 1) complex 
(nX 1) complex 
n v X (jl(j) + 77^) complex, eq(32) [6] 
n v X (n^ + 77^) complex, eq(33) [6] 
(rx 1) complex/real 
real scalar 
(r X 2) real matrix 

K + n v) 


24 


3.3.3 form sample 

form_sample 


Purpose 

Divide time record into windowed-over lapped-zero appended samples. 


Synopsis 

[ymeas_samp,yDFT_samp,omeg,Ts,t_samp,nout, nsamp] = f orm_sample (ymeas , t , . . . 
winpara, wintype) 


Description 

Given a discrete time record, ymeas as a 2-D matrix (time varying row-wise, different channel 
column- wise) and a column vector of discrete time, t, this routine construct samples with 
winpara(l) non-zero point lengths. The output variable, ymeas_samp, is a 3-D matrix (varying 
time, output #, sample #). The corresponding output variable, yDFT_samp denotes the DFT 
of ymeas_samp and stored as a (1 x nsamp) array of FRD models. The variable omeg denotes 
a column vector of frequencies corresponding to the DFT of the sample time vector, t_samp 
having a sampling interval Ts for nout channels. The integer nsamp denotes the number of 
samples generated. 

The input parameter winpara(2) specifies the total number of data points in each sample 
(which the FFT is applied). Zeros are appended if winpara( 2) > winpara( 1). The number 
of samples can be increased by specifying the number of overlap points, winpara(3). 

Each sample can be windowed by specifying the parameter, wintype as 1 rec 1 for rectangular, 
‘han’ for Hanning, and ‘ham’ for Hamming windows. The windowing is applied only to the 
nonzero data component of length winpara(l). 
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3.3.4 fungrad 

fungrad 


Purpose 

Defines the cost function and its gradient with respect to variables in a smallest model 
validating set optimization problem. 

Synopsis 

[f , g] = fungrad (x , xio , et ao , EF , PEFGH , Wunc , bo , blkf , cdimM , rankN , rankM , nrrs , nr cs ) 

Description 


This subroutine computes the cost function, f and its gradient, g, for a smallest set opti- 
mization problem which is used by the optimizer fmincon in MATLAB Optimization Toolbox. 
The above cost function and its gradient are evaluated at the given input optimization variable: 


x(l : n$) 



= Rety) 

( n ^ 

X 

1) 

real 

x(n </, + 1 : 



= Im{ij}) 

(n^ 

X 

1) 

real 

x(2n $ + 1 

• 277 

+ n^) 

= Re(<t>) 

( n <t> 

X 

1) 

real 

x(2 + n^ + l 

m) 

= Im{4>) 

( n <t> 

X 

1) 

real 

x(m + 1 : 

m + n 

rrs) 

= d \ , . . . , d nrra 

(fires 

X 

1) 

real 

x(m + n rrs + 1 : 

m + n rrs + 2 n rcs ) 

= Re(d nrrs+ 1 ), Im(d nrrs+1 ), . . . 








Rc{dn rrs +n rcs ) } I m {d nrrs +n r C3 ) 

(2n rcs 

X 

1) 

real 

x(m + n rrs + 2 n 

res R 1) 

= x 2 




real 


where 771 := 2 + 2 770. The remaining input arguments are subroutine parameters defined 
in the command description congrad. 
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3.3.5 mnmvcl 


mnmvcl 


Purpose 

This program sets up and solves for a smallest scaled model validating set based on [6]. 

Synopsis 

[unc , Jmin, options , vFLAG] = mnmvcl (rspec ,yspec ,Gaug,blkf , Vnoise , Vdist ,Ts , . . . 
Wunc , iv_low , i v Jiig , nskip , ndist ) 


Description 

Given DFTs of the command or plant (for open loop system) input, rspec and plant output, 
yspec, and allowances for the unknown exogenous output noise and disturbances defined by 
filters Vnoise and Vdist respectively, this program attempts to solve for a smallest scaled 
uncertainty norm bound unc amongst all feasible LFT plant sets satisfying model validation 
conditions in the frequency domain. The minimization is over all feasible LFT augmented 
plant sets defined by a nominal 2 by 3 block canonical augmented plant, Gaug, having the 
block uncertainty structure blkf . The augmented plants are assumed to be discrete linear time- 
invariant systems with a sampling interval Ts, identical to the original discrete time signal which 
produced rspec and yspec. 

Figure 7 shows the LFT augmented plant with all relevant signals required in the optimiza- 
tion problem. The block diagonal uncertainties are assumed to be in the forms: 

1! • diag(wi/ ni , • • ■ , w r I nr , w r -\-iIn r +i , • ,w r _^_ c I nr + ci x ' w r -\- c -\-iI nr ^ c ^ i: • • • ^X’W T I nT ^j 

^ B • — diag((5i/ ni , ■ ■ • , &rln r ? &r+lln r +i , , ^r-\-c^n r ^. c > ^r+c+ 1 , ' ' ; A r ) 

where Si E 71, for i = 1, . . . , r; Si G C, for i — r + 1, . . . , r + C; A * G C n ' lXm ‘, 
for i = r + C + 1, . . . , T. The symbol, A b, denotes a block diagonal uncertainty structure 
bounded by a unit ball. The symbols r, C, and T denote the number of repeated real scalar 
blocks, repeated complex scalar blocks, and total number of uncertainty blocks, respectively. In 
general, the scalar Si will be a function of frequency, in contrast to “parametric” uncertainties 
which are sometimes explicitly assumed to be constants. 

The frequencies, iv_low and iv_hig define the bandwidth over which model validation con- 
ditions are to be satisfied, nskip denotes the number of frequency points to be skipped when 
the frequency points based on fixed time sample period are overly dense, ndist denotes the 
number of unknown exogenous disturbance input. 

mnmvcl solves for a smallest scale model validating set based on the following nonlinearly 
constrained optimization formulation: 


subject to 


mm x~ 

4>£C n <i> rtec 71 * ,&i,...,5 r +c,x 2 en+ 


&(</>, i>) - Sii)^, VO = o, Si e Hi 


\Si 


IwA < 0 
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i — r + 1, . . -i , r + c 


6(0, 0) - A??i(0, 0) = 0, 6 e C,: 

|5j| 2 - |^i| 2 < o 

116(0, 0)|| 2 -x 2 |^| 2 ||r ?l (0,0)|| 2 < 0, i = r + c+1,...,t 

II0II 2 < bl 


Each element in Wp, . . . , W T of the user specified frequency varying column vector, Wunc, defines 
a radius for each individual uncertainty block. The scaling factor, X, appears only with the 
radii of the nonrepeated complex blocks, X ■ Wi,i = r + C + 1, . . . , T, implying that the 
scaling minimization is applied only to the nonrepeated complex blocks. This is an assumption 
in the problem statement which can be modified if additional a priori assumptions on the 
uncertainty structure indicates otherwise. We seek a smallest scaling, X, such that a set of unit 
norm-bounded model errors, <5i, . . . , 6+ci and A. r _|_ c+ i, . . . , A r (parameterized by 6 1/), will 
exactly reproduce the given input and output signals to within unknown but bounded exogenous 
signal allowance limits (see routine noise allow for a description of signal allowances). It is 
important to note that the smallest scaling is with respect to user defined radii for each individual 
nonrepeated complex uncertainty block. 

The output unc is a frequency varying column vector denoting a set of smallest model 
validating uncertainty norms satisfying 


unCi : = 



|6| < \wi\, i — 1, — ,r + c 


for the repeated real and complex blocks and 


urtCi 



1 1 Ai| | < x\wi\, 


i — r + c+1,. . . ,r 


The second output variable Jmin denotes the minimum cost function. Notice that satisfying the 
model validation conditons imply the existence of a set of 6 and Aj about a nominal model but 
in no way imply that the “true” model error has been recovered or that a “true” model even 
exists. 

mnmvcl performs a sufficiency test for convexity based on a norm test of the augmented plant 
blocks [Gil, G12] • If this test is passed the subroutine ubid_lmi is called which sets up a 
convex program and solves it using the subroutine mincx in LMI Toolbox. If the sufficient 
test fails, then mnmvcl assumes a general nonlinearly constrained optimization problem and 
uses a sequential quadratic programming algorithm as implemented in fmincon, a routine in 
OPTool. fmincon in turn calls fungrad which defines the cost function and its gradient, 
and congrad which defines the constraints and its gradients. 
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3.3.6 noise allow 


noise_allow 


Purpose 

This routine help generate noise (and disturbance) weight filters to unknown but bounded 
exogenous signal allowances used in model validation. 

Synopsis 

[vVnoise,Vnoise] = noise_allow(N,Ts,nout, label, ynoise,t_raw) 

Description 

The input variable N denotes an N-point DFT sample used in model validation with sampling 
interval Ts while the variable nout denotes the number of channels in the noise/disturbance 
signals. If a sample discrete time history, ynoise, is not available, its mean square value (or 
variance) of the signal is assumed known. The variance for each of the channels are entered 
interactively when the routine is called with only 3 or 4 input arguments, label is an optional 
script label. If a sample of discrete time noise sequence ynoise (2-D matrix with time varying 
row-wise) is known, the routine is called with all 6 input arguments including a column vector 
of discrete time points, t_raw. 

The output variable vVnoise denotes the frequency-varying noise weight in FRD format, 
computed from the fitted continuous state space filter Vnoise in SS format. 


Approach 

Given only the variance of a broadband (power content up to Nyquist frequency) discrete 
time random signal, we approximate a constant power spectra, S a from the relation [16] 

£ h 2 ] » 

-L s 

The above approximately constant power spectra can be related to a corresponding constant 
DFT magnitude using the periodogram formula 


So = S k = - |Xf| 2 , k=l,...,N 

where Xj} 7 denotes the DFT of an N-point sequence Uq, , V]\ T - \ 


N-l 

N \ - jlTTik 

x k : = Uie N 
i=0 

In summary, for the case when only the variances of an independent vector noise signal is known, 
its noise allowance in terms of its N-point DFT is defined by the constant matrix 


Vnoise 


NT s n ou t 

27rdiag (S[y 2 \ u . . . ,S[u 2 ] nout ,) 
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The factor \J'n 0U t is due to the unity 2-norm bound on the magnitude DFT of the vector signal 
at each frequency. This means that for multiple channels, a generous level of noise allowance 
can arise. 

For the case when a sample discrete time history ynoise which is significantly longer than 
N is known, several N-point subsamples can be partitioned and weighted (if desired) and an 
allowance is estimated by directly fitting a transfer function over the magnitudes of N-point 
DFT subsamples. The magnitude of DFTs are fitted interactively using stable real rational 
transfer function using the /i-Tool command drawmag. The input-output relationship is 

X* = V(u k )W?, k = 1, N/2 

where W f F denotes the N-point DFT of a signal Wq, ...,Wn- ± whose DFT is unity 2- norm 
bounded and V (cJk) denotes the frequency response of a noise allowance filter at frequency 
CPfc which maps a unity 2-norm bounded DFT signal to the given time sequence. Therefore, 
assuming a simulated noise whose DFT satisfies \Wf^ — 1 for all k, fitting | Xj ? | is equivalent 
to fitting \V (oUk)\ as required. 
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3.3.7 ovbndunc 


ovbndunc 


Purpose 

Overbound frequency varying uncertainty norm data using stable rational transfer functions. 


Synopsis 

[Wid_z] = ovbndunc (unc_frd,blkf , omeg, iv_low, iv_hig,Ts) 


Description 

The components of the FRD model unc_frd consists of the computed model validating 
uncertainty norms for each block whose structure is defined by each row of blkf . The input 
omeg is the linearly spaced frequency points from the original FFT of the discrete time data 
whose sampling interval is Ts seconds. These frequency points are used only to compare the 
rational function fit, Wid_z, to the data itself whose range is limited to the frequency bandwidth 
in radians defined by bounds iv_low and iv_high. The output, Wid_z, is a fitted discrete 
time state space model in SS format. The actual fitting is done interactively using the /i-Tool 
command, drawmag. 

For a tighter overbound, a linear programming approach given in [17] is recommended. 


31 


3.3.8 pinv frd 

pinv_frd 


Purpose 

Compute the pseudoinverse of a FRD model. 


Synopsis 

[piM_frd] = pinv_frd(M_frd,tol) 


Description 

This routine computes the pseudoinverse (or inverse), piM_frd of the input FRD model, 
M_frd. At each frequency point, the MATLAB function pinv is used with zero threshold 
defined by tol. 
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3.3.9 sseigunc 


sseigunc 


Purpose 

Form augmented state-space model with eigenvalue uncertainties. 

Synopsis 

[Prr ,blk_eig,nxi_eig,eigz, nreal] = sseigunc (Gnomd.Ts) 


Description 

Input Gnomd denote state space matrices for a discrete time nominal system with sampling 
interval Ts. The given nominal state matrix is block diagonalized to 2 by 2 blocks for complex 
conjugate pairs and 1 by 1 diagonal element for a real eigenvalue using blkd2x2 routine which 
also determines the number of purely real eigenvalues, nreal, placed first on the diagonal. 

When the subroutine is called the user is prompted for the set of nominal eigenvalues to 
be included as parametric uncertainties. The uncertainty for a purely real nominal eigenvalue 
is assumed to be purely real also while the uncertainties for each complex conjugate pair is 
also assumed to be complex conjugate, for consistency. All uncertainties in the eigenvalues 
are assumed to be independent and real and consists of the real perturbations in the purely 
real eigenvalue and the real and imaginary component perturbations in the complex conjugate 
eigenvalues. 

Prr denotes the resulting augmented plant in discrete time state space, SS format, as shown 
in figure below. Same notations are used as in fi - Tools User’s Guide [3] on pages 4-20 which 
describes a general affine perturbation for a linear state space model. The uncertainty structure 
definition matrix blk_eig reflects a set of real uncertainty variables which also includes repeated 
scalars if complex conjugate eigenvalues are present. The remaining output arguments include 
the number of channels due to eigenvalue uncertainties (i.e. dimension of T/ or £), nxi_eig, 
nominal eigenvalues in z-plane, eigz, and the number of purely real eigenvalues, nreal. 



% 

u 
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3.3.10 svd frd 


svcLfrd 


Purpose 

Compute the SVD of a FRD model and its rank at each frequency. 


Synopsis 

[uu,ss,vv,rankss] = svd_frd(M_frd,tol) 


Description 

This routine computes the SVD of the input FRD model, M_frd where uu*ss*vv = M_frd. 
At each frequency point, the MATLAB function svd is used with the rank threshold defined 
by tol. The output variable, rankss is a column vector of rank at each frequency. 
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3.3.11 ubid 


ubid 


Purpose 

Main program for computing a smallest model validating set. 


Synopsis 

This is a sample main program which calls the following subroutines to compute a smallest 
model validating set: 

form sample 
noise_allow 
ubid_uncstr 
weightunc 
mnmvcl 
ovbndunc 

Description 

To compute a smallest model validating set, various assumptions are made on the time 
signals to be used in defining the set of all model validating plants, the nominal model, the 
uncertainty structure, and exogenous noises and disturbances. These assumptions include: 

• N-periodic input/output discrete time signal samples 

• discrete linear time-invariant systems 

• unknown but bounded exogenous noise and disturbances 

• Linear Fractional Transformation (LFT) uncertainty structure 

• known relative importance of the initial uncertainty weights 


preprocess sample data 

define noise and disturbance allowances 

help construct uncertainty structure 

define relative uncertainty weight for optimization 

check feasibility and optimizes for a smallest set 

overbound uncertainty norms using stable rational filters 
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3.3.12 ubid lmi 


ubid_lmi 


Purpose 

Sets up convex program of a smallest set optimization problem. 

Synopsis 

[xopt , copt , EXITFLAG] = ubicLlmi (n_phi , n_psi , xio , et ao , coXI , Wunc , bo , blkf ) 


Description 


Sets up and solves a class of convex optimization problems defined by [Gn, G12] — 0 . For 
more general convex problems, this routine can be modified as necessary. The convex optimizer 
mincx found in LMI Toolbox [20] is used which assumes a linear cost function 

min c T z 

Z 

subject to LMI constraints 

%z 


where z := [Re(0); Im^); Re(0); lm(0); X 2 ] G ‘Jl 2n -<l>+ 2n <l>+ 1 _ The constants, C, Qi, Si, L, 
and A, are defined in equation 22. 

The input parameters, n_phi and n_psi denote the dimensions of the complex vectors, 0, 
and ip, which are the free parameters defining all model validating sets, xio and etao denote 
the constant values of xi and eta which are the fictitious complex vector signals going out and 
into the uncertainty block, respectively. The complex matrix, coXI denotes the coefficient of 
the vector (0; 0) in ^-equation, i.e. , [Oj ; . . . ; f2 r ] where Qi is defined in equation 13. The 
column vector Wunc denotes the uncertainty weights used in the search for a smallest model 
validating set. It is used to define the relative significance of the various uncertainties which 
together define a particular model validating set. bo denotes the exogenous signal uncertainty 
allowance radius, blkf defines the uncertainty structure. 

The output variable x = xopt are defined by: 


- H6vll 2 

SiZ 


sym 

I 


b 2 0 syn 

Lz I 


> 0, i — 1, . . . , r 

> 0 


x(l : n_psi) 

x(n_psi+l : 2n_psi) 

x(2n_psi+l : 2n_psi+n_phi) 

x(2n_psi+n_phi+l : 2n_psi+2n_phi) 

x(2n_psi+2n_phi+l) 

copt 

EXITFLAG 


= Re(0) 

= Im (0) 

= Re(0) 

= lm(0) 

= X 2 

= global minimum 
= 1, global minimum 

0, convex program infeasible 
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3.3.13 ubid uncstr 


ubicLuncstr 


Purpose 

Form augmented plant with uncertainty structure. 


Synopsis 

[Paug, blkf ,dimp,uncstr_ch.oice] = ubid_uncstr(Po,Ts) 


Description 


Given a nominal model, Po this routine help to construct an augmented plant, Paug with 
the uncertainty structure, blkf in the following canonical form: 



Gil 

G21 


G 12 

G 22 


G 13 
G23 



where £, = A 77. The specific forms for G and A depends on the particular type of uncertainty 
structure selected, uncstr_choice, from the following short list of commonly used structures: 


uncstr_choice 

Type of uncertainty 

0 

provided by user (filename, all output variables defined in Synopsis) 

1 

additive 

2 

additive + eigenvalue 

3 

full output multiplicative 

4 

full output multiplicative + eigenvalue 

5 

diagonal output multiplicative 

6 

diagonal output multiplicative + eigenvalue 

7 

full input multiplicative 

8 

full input multiplicative + eigenvalue 

9 

diagonal input multiplicative 

10 

diagonal input multiplicative + eigenvalue 


The six integers, dimpil) . . . , dimp{ 6), define the dimensions of the corresponding vectors, 
V, £, V’ A A and r. 

Alternately, the user can load the augmented plant and uncertainty structure from a file by 
choosing uncstr choice = 0 option. In this case all variables, Paug, blkf, and dimp must be 
provided. 

The input variable Ts denotes the sampling interval in seconds and is necessary only if 
eigenvalue uncertainties are included. 
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3.3.14 weightunc 

weight unc 


Purpose 

Define initial uncertainty weight/radii for smallest set optimization. 


Synopsis 

[Wunc ,par_bnd] = weightunc (blkf ,omeg,Ts) 


Description 

This routine helps the user to select uncertainty weights for use during optimization for 
finding a smallest scaled model validating set. The selection of these weights should be based on 
the relative importance of different uncertainty blocks and over frequency for a given uncertainty 
block. This weighting reflects a priori knowledge by the user on the uncertainty structure of the 
system, however, the weighting over frequency has not been implemented yet. The optimization 
problem involves finding a minimum scaling of these user defined weights and still satisfy model 
validation conditions. 

The input consists of the uncertainty block structure, blkf using /i-Tool convention [3], the 
discrete frequency points where weighting needs to be defined, omeg, and the sampling time of 
data, Ts. The output Wunc, is a frequency varying column vector uncertainty weight in FRD 
format, to be used in a numerical optimization for a smallest set while par_bnd defines a column 
vector of parameter uncertainty allowance. 

A basic question of whether a particular set of plant model defined by Wunc cannot be 
invalidated by a given set of data is a special case of this smallest scaled set optimization 
problem. Namely, if Jmin < 1 (from subroutine mnmvcl) then the given set of plant model 
defined by Wunc cannot be invalidated by the given data. If the uncertainty structure leads 
to a convex program (including a special class where \_G\i-,G\ 2 \ — 0) then, Jmin < 1 is a 
necessary and sufficient test. For a more general uncertainty structure where Jmin is the result 
of a nonliearly constrained optimization problem, the above unit inequality becomes effectively 
a sufficient test since it is usually difficult to guarantee that a global optimum has been attained. 
Hence if Jmin exceeds 1, the test result is indeterminate. 
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4 Application Example 1: Uncertainty bounds of a large flexible 
structure 

In this application example, we consider the problem of estimating model uncertainty bounds 
from data for active vibration control of a large flexible structure. We also show in detail 
how UBID Toolbox can be used in this problem. 


4.1 Control-Structure-Interaction Evolutionary Model 

Figure 11 is a schematic of the NASA Langley’s flexible structure testbed. We consider a 
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Figure 11: Control-Structure-Interaction Evolutionary Model (CEM). 


NASTRAN model of the CEM consisting of 638 finite element grid points and include only 
the first 30 structural modes which are below 8 Hertz. The expected controller bandwidth is 
about 2 Hertz and the structural responses are sampled at 20 Hertz. Although many possible 
actuator and sensors are available for control as shown in the figure, we consider only 3 air 
thrusters and 3 displacement sensors which are approximately collocated and colinear in the 
x, y, and z axis corresponding to directions 5, 1, and 2 in the Figure 11. For more details 
please refer to [18]. 

For the purposes of illustration and evaluation of the uncertainty bound estimation ap- 
proach, we assume that a physical system and associated exogenous noises can be represented 
by a truth model consisting of a 30 structural mode NASTRAN model. We denote the 60-th 
order true state space model discretized at sampling interval T s as 

A-true B f %k\ 

c D \ \u k ) 

Response measurements are simulated using this truth model and are used in estimating the 
uncertainty bounds. 
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4.2 Nominal model and uncertainty structure 


Consider a perturbed system is defined by ( A pert , B, C, D ) where 

■Apert ■ A-true "F A A 

such that Atme and A A matrices have the following block diagonal forms 



' Ai 

0 1 




■A-true • 

0 

1 

o 

CO 

<1 

e ^ 60x6 °, 

Ai := 




’ A A, 

0 ' 





AA := 

£A 3 o . 

e ft 60x60 , 

SA t := 

5ai 

. SPi 

<5 A ' 

5oti 


0 






i — 1, . . . , 30 


i = 1, . . . , 30 


Denoting the true eigenvalue conjugate pairs as A* ± := ± /3*. its perturbations denoted 

as 5Aj ± are assumed constants such that A pert is stable, i.e., |Aj ± | < 1, |Aj ± + <5Aj ± | < 1, 
i = 1, . . . , 30. In this example, all eigenvalues are assumed to occur in complex conjugate 
pairs for simplicity although in general eigenvalues on the real axis in the z-plane can also 
be included. 


To generate a nominal model as typical in a realistic application, we internally balance 
the perturbed system containing erroneous eigenvalues, (A pert , B,C, D), using a similarity 
transformation, P , so that 


P^ 1 A pcrt P P~ X B 1 (X k \ 

CP D \ [uj 

where x k Px k denotes the balanced state at time k. The nominal model is then obtained 
by keeping only the states corresponding to the 12 largest Hankel singular values. The 
neglected 48 states reflect unmodeled dynamics. 

The above procedure to generate a nominal model is intended to reflect a common ap- 
plications scenario, where various finite element models with a large number of states are 
available for a structural system, whos eigenvalues are not totally reliable. The nominal 
plant model is derived from the perturbed 60 state model, whos frequencies and damping 
values differ from the truth model even before an order reduction. This reflects inevitable 
parametric uncertainties about nominal eigenvalues. In addition, notice that the true mode 
shapes are preserved in the reduced nominal model to the extent of the approximation of the 
principal component vectors to the modal coordinates for lightly damped flexible structures 
19]. Interestingly, in addition to the fact that eigenvalues define the frequency and damping 
of resonances in structural vibrations, eigenvalues as candidates for parametric uncertain- 
ties for flexible structures appears huristically attractive because they are also invariant to 
similarity transformation. 

For the above reasons, we select a subset of the eigenvalues from the nominal model 
as candidates for parametric uncertainties and additive uncertainties to take into account 
unmodeled dynamics and mode shape errors, as illustrated in Figure 12. To demonstrate 
that this choice in the uncertainty structure is sufficient, suppose a true system consists of 
a high order state space system (A, B,C, D) where A has 2 by 2 block diagonal structure, 
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Figure 12: Uncertainty structure connections for a flexible structure 
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without loss of generality. Furthermore, suppose the matrices are partitioned as follows: 



A\ 



' Ai + 8 A x 


A2 



A 2 + 8A 2 

B i ' 


' B 1 

+ 8B1 ' 


B 2 . 


b 2 

T 8B 2 _ 


C x 

C 2 ] 

■■=[ 

C x + 8Cx C 2 + 5C 2 } 


D := D + 8D 

where (Ai, D) is intended to represent a lower order nominal model. Then the in- 

put/output mapping of the true system can be written as 


where 


/ 


nominal model 
with eigenvalue error 

A 


\ 


y — C\(zl — Ai — 8Ai) Bi + D +A add 


u 


V 


A 


add 



f eedthrough 
gain error 


modal output error unmodeled dynamics 

.A. ✓Si 

+ SCi(zI — Ai) + Ci(zl — A\) ^6Bi -\- C2{zl — A2) 1 B2 

^ V 

modal input error 


The point here is that such true model can be captured by additive uncertainty structure, 
A add , about a nominal model with eigenvalue uncertainties. The subsequent issue is how 
large are the uncertainty balls around the nominal model. Note that the modal output and 
input error terms arise from the errors 8C\ and 8B\ about their respective nominal values 
C\ and B\. Alternatively, 8C\ and 8B 1 (along with 8D) could be viewed and treated as 
additional parametric uncertainties in a state space model similar to the treatment of eigen- 
value uncertainties captured in 8A\. However, the former untertainty terms are coordinate 
dependent unlike eigenvalues and hence difficult to treat them as parametric uncertainties. 

Figure 13 shows the eigenvalues of the true and nominal models. The bottom figure is 
a magnified view of the 3 dominating structural resonance modes. The true and nominal 
eigenvalues appears to match quite closely although the subsequent frequency response com- 
parisons show significant differences. The above tasks to generate the true and nominal 
models are given in subroutine ex_cem_data. 

A comparison of the magnitude gains between the true (solid) and nominal (dash) models 
and the corresponding differences (dot) for all inputs and outputs are shown in Figure 14. The 
difference plot (dot) show significant discrepencies especially around structural resonance 
frequences. For example, consider a Bode plot for input 1 to output 1 as given in Figure 
15. The figure shows that despite significant errors near their resonance frequencies, the 
nominal model approximates the three low frequency structural resonances at .9, 4.5 and 
15 rad/sec Several higher frequency secondary resonances beyond 20 rad/sec are ignored by 
the nominal model. The largest discrepency between the true and nominal model frequency 
responses occurs at the dominating structural resonance at .9 rad/sec. 


4.3 Smallest unmodeled dynamics for model validation 

The main goal is to determine if the uncertainty bound estimation algorithm can capture both 
the true unmodeled dynamics and parametric uncertainties from data with reasonable a priori 
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CEM eigenvalues: True (o) vs Nominal (x) 




Re(z) 


Figure 13: True and nominal eigenvalues of CEM model. 


assumptions including noise allowances to account for unknown exogenous disturbances. In 
the following, we show the results based on the simplest case where only additive uncertainty 
is used to two other cases where we incorporate a priori parametric eigenvalue uncertainty 
allowances while seeking smallest additive uncertainties for model validation. 

All three input test signals to the true system is generated by independently driving a 
fourth order Butterworth filter of bandwidth 5 Hertz with a white (bandwidth of Nyquist 
frequency) Gaussian random signal with unit standard deviation. The additive output noise 
is similarly chosen to be a white Gaussian random signal with a standard deviation of .01. 
The test data consists of two periods each about 408 seconds duration. The test signal is 
only applied during the first 204 seconds followed by free response for another 204 seconds. 
Notice that the next period begins after most of the transient free response decays to the 
output noise levels. This is designed to mitigate the erronous effects of nonperiodicity in the 
computation of DFT’s for use in model validation. 
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Output 3, Input 1 Output 2, Input 1 Output 1 , Input 1 








Figure 14: Frequency response magnitude compari 


, Difference(dot) 



rad/sec 


of CEM: outputs 1-3, inputs 1-3 












4.3.1 Smallest additive uncertainty 


Step 1 Load output and input time histories 

• Load raw time histories ymeasuraw, rmeas_raw: 

LOAD SIMULATED INPUT/OUTPUT DATA, ymeas_raw, rmeas_raw, t_raw 
>> load ex_cem_data.mat 




Time (sec) 


• Window data to form multiple samples by specifying number of nonzero points in 
sample, npts, number of points in sample to apply FFT, nfft, number of data 
points to overlap, noverlp. Plots for the sample time histories to be used in 
uncertainty bound estimation 

>> [ymeas,yDFT,omeg,Ts,t,nout,nsamp] = f orm_sample(ymeas_raw,t_raw,winpara, "rec") 
PREPARE TIME RECORDS (WINDOWED-OVERLAPPED-ZERO APPENDED SAMPLES) 

Beginning time of data available = 0 

End time of data available = 819.15 

Number of time points in signal sample (npts) = 8192 

Number of points in sample (nfft), with zeros appended = 8192 

Number of data overlap (noverlp) = 0 

>> [rmeas,rDFT,omeg,Ts,t,ninp,nsamp] = f orm_sample(rmeas_raw,t_raw,winpara, "rec") 
PREPARE TIME RECORDS (WINDOWED-OVERLAPPED-ZERO APPENDED SAMPLES) 

Beginning time of data available = 0 

End time of data available = 819.15 

Number of time points in signal sample (npts) = 8192 

Number of points in sample (nfft), with zeros appended = 8192 

Number of data overlap (noverlp) = 0 

There are 2 data samples - select sample number for MV.. 1 
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CSI EVOLUTIONARY MODEL data, Ts=0.05, Sample 1, 8192 points 



Time (sec) 



Step 2 Form allowance filters for output noise and exogenous disturbances 

• Load output noise time history variable, ynoise and fit filter, Vnoise 

LOAD TRUE OUTPUT NOISE, ynoise 
>> load ex_cem_truenoise .mat 

>> [vVnoise , Vnoise] = noise_allow(ntime ,Ts ,nout , ["noise_allow_" label ynoise, t_raw) 

PREPARE TIME RECORDS (WINDOWED-OVERLAPPED-ZERO APPENDED SAMPLES) 

Beginning time of data available = 0 

End time of data available = 819.15 

Number of time points in signal sample (npts) = 8192 

Number of points in sample (nfft) , with zeros appended = 8192 

Number of data overlap (noverlp) = 4096 

Interactively fit magnitude DFT, Channel 1.. 

Re-fit Channel 1 ? (y/n) n 
Interactively fit magnitude DFT, Channel 2.. 

Re-fit Channel 2 ? (y/n) n 
Interactively fit magnitude DFT, Channel 3.. 

Re-fit Channel 3 ? (y/n) n 
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• Define allowance filter for exogenous disturbances, Vdist from expected variance 

>> [vVdist , Vdist] = noise_allow(ntime ,Ts ,ninp, ["dist_allow_" label 
CONSTRUCT EXOGENOUS SIGNAL ALLOWANCE FILTER FROM VARIANCE ESTIMATE 

Enter variance (MSV) of unknown signal, Channel 1.. 0 

Enter variance (MSV) of unknown signal, Channel 2.. 0 

Enter variance (MSV) of unknown signal, Channel 3.. 0 


Step 3 Load nominal model, Po 

LOAD NOMINAL MODEL, Po (SS format) 

>> load ex_cem_nominal.mat 

TRUE EIGENVALUE PERTURBED EIGENVALUE DIFFERENCE 

ans = 


0.9954 

+ 

0 . 0462i 

0.9943 

+ 

0 . 0476i 

-0.0011 

-0.0015 

0.9951 

+ 

0 . 0467i 

0.9953 

+ 

0 . 0470i 

0.0001 

-0.0003 

0.9940 

+ 

0 . 0482i 

0.9948 

+ 

0 . 0474i 

0.0009 

0.0008 

0.9689 

+ 

0 . 2280i 

0.9703 

+ 

0 . 2264i 

0.0014 

0.0016 

0.9673 

+ 

0 . 2342i 

0.9673 

+ 

0 . 2337i 

-0.0000 

0.0006 

0.9555 

+ 

0 . 2754i 

0.9560 

+ 

0 . 2751i 

0.0005 

0.0003 

0.8830 

+ 

0 . 4490i 

0.8827 

+ 

0 . 4485i 

-0.0003 

0.0005 

0.7363 

+ 

0 . 6552i 

0.7364 

+ 

0 . 6549i 

0.0001 

0.0003 
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0.7109 

+ 

0 . 6816i 

0.7114 

+ 

•H 

CT> 

O 

00 

QO 

O 

0.0005 

0.0006 

0.7021 

+ 

0 . 6903i 

0.7019 

+ 

0 . 6882i 

-0.0001 

0.0021 

0.6758 

+ 

0 . 7150i 

0.6766 

+ 

0.7148i 

0.0008 

0.0002 

0.5948 

+ 

0 . 7718i 

0.5962 

+ 

0.7733i 

0.0014 

-0.0014 

0.5179 

+ 

0 . 8319i 

0.5180 

+ 

0 . 8324i 

0.0002 

-0.0005 

0.5155 

+ 

0 . 8333i 

0.5148 

+ 

0 . 8351i 

-0.0007 

-0.0018 

0.3825 

+ 

0 . 8989i 

0.3822 

+ 

0 . 9003i 

-0.0003 

-0.0014 

0.2570 

+ 

0 . 9397i 

0.2590 

+ 

0 . 9392i 

0.0020 

0.0005 

0.2040 

+ 

0 . 9515i 

0.2042 

+ 

0 . 9522i 

0.0002 

-0.0006 

0.0021 

+ 

0 . 9691i 

0.0013 

+ 

0 . 9675i 

-0.0008 

0.0017 

-0.0630 

+ 

0 . 9658i 

-0.0637 

+ 

0 . 9658i 

-0.0007 

-0.0000 

-0.0763 

+ 

0 . 9645i 

-0.0758 

+ 

0 . 9666i 

0.0005 

-0.0021 

-0.0837 

+ 

0 . 9638i 

-0.0838 

+ 

0 . 9636i 

-0.0001 

0.0001 

-0.0971 

+ 

0 . 9622i 

-0.0965 

+ 

0 . 9602i 

0.0006 

0.0020 

-0.1948 

+ 

0 . 9453i 

-0.1939 

+ 

0 . 9450i 

0.0009 

0.0003 

-0.2219 

+ 

0 . 9387i 

-0.2208 

+ 

0 . 9394i 

0.0011 

-0.0007 

-0.3244 

+ 

0 . 9061i 

-0.3240 

+ 

0 . 9049i 

0.0004 

0.0013 

-0.3276 

+ 

0 . 9049i 

-0.3292 

+ 

0 . 9039i 

-0.0016 

0.0010 

-0.3396 

+ 

0 . 9002i 

-0.3394 

+ 

0 . 8991i 

0.0002 

0.0010 

-0.4943 

+ 

0 . 8214i 

-0.4933 

+ 

0 . 8220i 

0.0010 

-0.0006 

-0.5133 

+ 

0 . 8091i 

-0.5131 

+ 

0 . 8083i 

0.0002 

0.0008 

-0.7389 

+ 

0 . 6002i 

-0.7396 

+ 

0 . 6018i 

-0.0007 

-0.0015 


NOMINAL EIGENVALUE (BALANCED REDUCED PERTURBED) 


ans = 


0.9950 + 0 . 0472i 
0.9950 - 0 . 0472i 
0.9943 + 0 . 0476i 
0.9943 - 0 . 0476i 
0.9701 + 0 . 2264i 
0.9701 - 0 . 2264i 
0.8838 + 0 . 4457i 
0.8838 - 0 . 4457i 
0.7121 + 0 . 6800i 
0.7121 - 0 . 6800i 
0.6759 + 0.7144i 
0.6759 - 0.7144i 
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lm(z) 


CSI EVOLUTIONARY MODEL eigenvalues: True (o) vs Nominal (x) 



Step 4 Choose additive plus eigenvalue uncertainty structure to form augmented plant 
Gaug. The first three modes are chosen as uncertain parameters. 

>> [Gaug.blkf ,dimp,uncstr_choice] = ubid_uncstr(Po,Ts) 

Form augmented plant with uncertainty structure 
OUTPUTS INPUTS 


eta 

(dimp(l) ) 

< 

[Gil , G12 , 

G13] 

< xi 

(dimp(2) ) 

y 

(dimp(3)) 

< 

[G21 , G22 , 

G23] 

< epsil 

(dimp(4) ) 






< nu 

(dimp(5) ) 






< u 

(dimp(6) ) 

ALL 

WEIGHTS FOR 

UNKNOWN 

EXOGENOUS 

INPUTS, Vdist=Vnoise 

! = I 


Select from the following list of uncertainty structures: 
AUGMENTED PLANT WITH . . . 

0. PROVIDED BY USER 

1. ADDITIVE UNCERTAINTY 

2. ADDITIVE + EIGENVALUE UNCERTAINTIES 

3. FULL OUTPUT MULTIPLICATIVE UNCERTAINTY 

4. FULL OUTPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

5. DIAGONAL OUTPUT MULTIPLIATIVE UNCERTAINTY 

6. DIAGONAL OUTPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

7. FULL INPUT MULTIPLICATIVE UNCERTAINTY 

8. FULL INPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

9. DIAGONAL INPUT MULTIPLIATIVE UNCERTAINTY 

10. DIAGONAL INPUT MULTIPLIATIVE + EIGENVALUE UNCERTAINTIES 
Enter the uncertainty structure desired (0 - 10) > 1 


ADDITIVE UNCERTAINTY 


eta_add(3) 

< 

— [Gil , G12 , G13] 

< 

xi_add(3) 

y(3) 

< 

— [G21 , G22 , G23] 

< 

epsil (3) 




< 

nu(3) 




< 

u(3) 
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Step 5 Choose reference weights, Wunc used in smallest set uncertainty optimization. As 
defined by blkf , the uncertainty blocks consists of six repeated real 2 by 2 blocks with 
a single full 3 by 3 complex block. 

>> [Wunc , par _bnd] = weightunc (blkf ,om,Ts) 


Define optimizing weights for uncertainty blocks . . . 
blkf = 


3 3 

Enter CONSTANT weights (over all frequencies) for each uncertainty block as a 
COLUMN vector (corresponds to BLKF variable) > 1 

Step 6 Solve for a smallest scaled model validating set. 

>> [unc_frd,x2, options, vFLAG] = mnmvcl (rDFTi ,yDFTi ,Gaug,blkf , vVnoise , vVdist ,Ts , . . . 
>> Wunc , iv_low , iv_hig , nskip , ndist ) 

Check feasible conditions (existence) for model validation: 

-Necessary condition eq(16) satisfied: M full row rank => eyo E Im(M) 
-Necessary condition eq(20) satisfied 
NM_beta is full row rank (T2 does not exist) => LHS of eq(20) is 0 


[Gil G12] = 0 satisfied, Convex Program solved using LMI Toolbox 


>> [x,copt,EXITFLAG] 
FREQ 1 0.107 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 2 0.199 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 3 0.291 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 4 0.383 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 5 0.476 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 6 0.568 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 7 0.660 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 8 0.752 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 9 0.844 (r/s) 

>> [x,copt,EXITFLAG] 
FREQ 10 0.936 (r/s) 


= ubid_lmi (n_phi , 
UNC = [ 0.002] 

= ubid_lmi (n_phi , 
UNC = [ 0.004] 

= ubid_lmi (n_phi , 
UNC = [ 0.003] 

= ubid_lmi (n_phi , 
UNC = [ 0.031] 

= ubid_lmi (n_phi , 
UNC = [ 0.021] 

= ubid_lmi (n_phi , 
UNC = [ 0.042] 

= ubid_lmi (n_phi , 
UNC = [ 0.025] 

= ubid_lmi (n_phi , 
UNC = [ 0.205] 

= ubid_lmi (n_phi , 
UNC = [ 0.419] 

= ubid_lmi (n_phi , 
UNC = [ 1.357] 


_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii ,etao_ii ,EF, 
_psi ,xio_ii ,etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 
_psi ,xio_ii , etao_ii ,EF, 


ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 
ii ,Wunc_ii ,bo ,blkf ) 


>> [x,copt,EXITFLAG] = ubid_lmi (n_phi ,n_psi , xio_ii , etao_ii , EF_ii , Wunc_ii ,bo ,blkf ) 
FREQ 200 18.423 (r/s) UNC = [ 0.000] 

>> [x,copt,EXITFLAG] = ubid_lmi (n_phi ,n_psi , xio_ii , etao_ii ,EF_ii , Wunc_ii ,bo ,blkf ) 
FREQ 201 18.515 (r/s) UNC = [ 0.000] 


51 


>> [x,copt,EXITFLAG] = ubid_lmi (n_phi ,n_psi , xio_ii , etao_ii , EF_ii , Wunc_ii ,bo ,blkf ) 
FREQ 202 18.607 (r/s) UNC = [ 0.000] 

>> [x,copt,EXITFLAG] = ubid_lmi (n_phi ,n_psi , xio_ii , etao_ii ,EF_ii , Wunc_ii ,bo ,blkf ) 
FREQ 203 18.699 (r/s) UNC = [ 0.000] 

>> [x,copt,EXITFLAG] = ubid_lmi (n_phi ,n_psi , xio_ii , etao_ii ,EF_ii , Wunc_ii ,bo ,blkf ) 
FREQ 204 18.791 (r/s) UNC = [ 0.000] 

>> [Wid_z] = ovbndunc(unc_frd,blkf ,omeg,iv_low,iv_hig,Ts) 

Curve fit uncertainty bound # 1, blkf = [3,3] .. 

Plot summary of weights? (y/n) n 


Figure 16 shows the smallest additive uncertainty (circle) that satisfy model validation 
conditions as compared to the true additive uncertainty (dotted line). A fitted additive 
uncertainty model is also shown (dashdot line) as a reference. The identified smallest additive 
uncertainty accurately recovers the true additive uncertainty. However, notice that the true 
additive uncertainty is due to a combination of both parametric eigenvalue uncertainties 
and unmodeled dynamics while the identified uncertainty model consists only of additive 
uncertainty. This demonstrates a limiting property that model validating conditions are 
only necessary conditions in determining “correct” uncertainty models, i.e. , just matching 
input-output responses are not germane to uncertainty structure. In this example, it is clear 
that a nominal allowance in the output noise did not have a significant effect on the identified 
additive uncertainty. Finally, note that due to the simplicity of the uncertainty structure, 


Case 1 



Figure 16: Identified smallest additive uncertainty (circle), true additive uncertainty (dot) 
the numerical computations involved convex optimizations and was very efficient. 
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4.3.2 Smallest additive uncertainty with eigenvalue uncertainty (.001) allowance 

Step 4 Choose additive plus eigenvalue uncertainty structure to form augmented plant 
Gaug. The first three modes are chosen as uncertain parameters. 


Form augmented plant with uncertainty structure 
OUTPUTS INPUTS 


eta 

(dimp(l) ) 

< 

[Gil, G12 , 

G13] 

< 

xi 

(dimp(2) ) 

y 

(dimp(3) ) 

< 

[G21 , G22 , 

G23] 

< 

epsil 

(dimp(4) ) 






c 

nu 

(dimp(5) ) 






< 

u 

(dimp(6) ) 

ALL 

WEIGHTS FOR 

UNKNOWN 

EXOGENOUS 

INPUTS, Vdist= 

=Vnoise 

— I 


Select from the following list of uncertainty structures: 
AUGMENTED PLANT WITH . . . 

0. PROVIDED BY USER 

1. ADDITIVE UNCERTAINTY 

2. ADDITIVE + EIGENVALUE UNCERTAINTIES 

3. FULL OUTPUT MULTIPLICATIVE UNCERTAINTY 

4. FULL OUTPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

5. DIAGONAL OUTPUT MULTIPLIATIVE UNCERTAINTY 

6. DIAGONAL OUTPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

7. FULL INPUT MULTIPLICATIVE UNCERTAINTY 

8. FULL INPUT MULTIPLICATIVE + EIGENVALUE UNCERTAINTIES 

9. DIAGONAL INPUT MULTIPLIATIVE UNCERTAINTY 

10. DIAGONAL INPUT MULTIPLIATIVE + EIGENVALUE UNCERTAINTIES 
Enter the uncertainty structure desired (0 - 10) > 2 


Mode number z-eigenvalue 


s-eigenvalue 


1 

2 

3 

4 

5 

6 


0. 99504+0. 047202i 
0. 99433+0. 047632i 
0. 97008+0. 22638i 
0. 88381+0. 44571i 
0. 71212+0. 68002i 
0. 67588+0. 71436i 


-0. 077057+0. 94803i 
-0. 090722+0. 95733i 
-0. 077366+4. 5851i 
-0. 2043+9. 3416i 
-0. 309306+15. 2468i 
-0. 33432+16. 2614i 


SELECT EIGENVALUES TO INCLUDE AS PARAMETRIC UNCERTAINTY 
Enter a row of Mode numbers (in increasing order) . . 


[1:3] 


AUGMENTED PLANT WITH EIGENVALUE UNCERTAINTIES 
State-space model with 15 outputs, 15 inputs, and 12 states. 
OUTPUTS 

due to 0 real eigenvalue uncertainties 

due to 3 complex conjugate eigenvalue pair uncertainties 
due to 3 physical outputs 
INPUTS 

due to 0 real eigenvalue uncertainties 

due to 3 complex conjugate eigenvalue pair uncertainties 
due to 3 physical inputs 

>> [Prr ,blk_eig,nxi_eig,eigz,nreal] = sseigunc(Po,Ts) 


AUGMENTED PLANT WITH ADDITIVE + EIGENVALUE UNCERTAINTIES 

eta_eig(12) < [Gil, G12, G13] < xi_eig(12) 

eta_add(3) < [G21, G22, G23] < xi_add(3) 
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y(3) 

< 

epsil (3) 


< 

nu(3) 


< 

u(3) 


Step 5 Choose reference weights, Wunc used in smallest set uncertainty optimization. As 
defined by blkf , the uncertainty blocks consists of six repeated real 2 by 2 blocks with 
a single full 3 by 3 complex block. 

>> [Wunc , par_bnd] = weightunc (blkf ,om,Ts) 

Define optimizing weights for uncertainty blocks . . . 

blkf = 

-2 0 

-2 0 

-2 0 

-2 0 

-2 0 

-2 0 

3 3 

Enter CONSTANT weights (over all frequencies) for each uncertainty block as a 
COLUMN vector (corresponds to BLKF variable) > [ . 001*ones (6 , 1) ; 1] 

Step 6 Solve for a smallest scaled model validating set. 

>> [unc_frd,x2, options, vFLAG] = mnmvcl (rDFTi ,yDFTi ,Gaug,blkf , vVnoise , vVdist ,Ts , . . . 

>> Wunc , iv_low, iv_hig,nskip ,ndist) 

Check feasible conditions (existence) for model validation: 

-Necessary condition eq(16) satisfied: M full row rank => eyo E Im(M) 

-Necessary condition eq(20) satisfied 
NM_beta is full row rank (T2 does not exist) => LHS of eq(20) is 0 


[Gil G12] .ne. 0, use SQP routine in Optimization Toolbox 


>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 1 0.107 (r/s) UNC = [ 0.003] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii , bo, blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 2 0.199 (r/s) UNC = [ 0.003] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 3 0.291 (r/s) UNC = [ 0.004] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii , bo, blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 4 0.383 (r/s) UNC = [ 0.008] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii , bo, blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 5 0.476 (r/s) UNC = [ 0.006] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii , bo, blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
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FREQ 6 0.568 (r/s) UNC = [ 0.002] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 7 0.660 (r/s) UNC = [ 0.002] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 8 0.752 (r/s) UNC = [ 0.010] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 9 0.844 (r/s) UNC = [ 0.059] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 10 0.936 (r/s) UNC = [ 0.154] MAX NUMBER ITERATION REACHED 


>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 200 18.423 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 201 18.515 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii , etao_ii ,EF_ii ,G0MEG_ii ,Wunc_ii , bo , blkf ,ncolM,minrankNMb ,minrankM,nrrs ,nrcs) 
FREQ 202 18.607 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 203 18.699 (r/s) UNC = [ 0.000] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[] ,vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 204 18.791 (r/s) UNC = [ 0.001] 

>> [Wid_z] = ovbndunc(unc_frd,blkf ,omeg,iv_low,iv_hig,Ts) 

Curve fit uncertainty bound # 1, blkf = [-2,0] .. 

Curve fit uncertainty bound # 2, blkf = [-2,0] .. 

Curve fit uncertainty bound # 3, blkf = [-2,0] .. 

Curve fit uncertainty bound # 4, blkf = [-2,0] .. 

Curve fit uncertainty bound # 5, blkf = [-2,0] .. 

Curve fit uncertainty bound # 6, blkf = [-2,0] .. 

Curve fit uncertainty bound # 7, blkf = [3,3] .. 

Plot summary of weights? (y/n) n 


Figure 17 shows the identified smallest additive uncertainty (circle) subject to eigenvalue 
uncertainty allowance (of .001 on the first three structural modes) as compared to the true 
additive error (dot). A fitted additive uncertainty model is also shown (dashdot line) as a 
reference. The level of the identified smallest additive uncertainty is smaller as expected 
than the true additive uncertainty because of the eigenvalue uncertainty allowance of .001 
on the first three structural modes. Apparently the eigenvalue uncertainty allowance of .001 
is not insignificant but less than the true eigenvalue error which is between .0003 to .0021 
(please refer to the eigenvalues given in Step 3 of the previous case). Obviously, there are 
no clear rules in the selection of a priori levels of parametric uncertainty allowances. How- 
ever, it is clear that incorporating “reasonable” levels of parametric uncertainty allowances 
may be helpful in determining accurate uncertainty models. In this case involving repeated 
scalar parametric uncertainties, solving a sequence of nonlinear optimization problems were 
required and the numerical computations involved were non trivial but is not central to the 
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Case 2 



Frequency (rad/sec) 


Figure 17: Identified smallest additive uncertainty (circle) with eigenvalue uncertainty allowance 
(.001), true additive uncertainty (dot) 

study reported here. 


4.3.3 Smallest additive uncertainty with eigenvalue uncertainty (.002) allowance 

Step 5 Choose reference weights, Wunc used in smallest set uncertainty optimization. As 
defined by blkf , the uncertainty blocks consists of six repeated real 2 by 2 blocks with 
a single full 3 by 3 complex block. 

>> [Wunc , par_bnd] = weightunc (blkf ,om,Ts) 

Define optimizing weights for uncertainty blocks . . . 

blkf = 

-2 0 

-2 0 

-2 0 

-2 0 

-2 0 

-2 0 

3 3 

Enter CONSTANT weights (over all frequencies) for each uncertainty block as a 
COLUMN vector (corresponds to BLKF variable) > [ . 002*ones (6 , 1) ; 1] 
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Step 6 Solve for a smallest scaled model validating set. 


>> [unc_frd,x2 , options , vFLAG] = mnmvcl (rDFTi ,yDFTi ,Gaug,blkf , vVnoise , vVdist ,Ts ,Wunc , iv_low, iv_hig, 
Check feasible conditions (existence) for model validation: 

-Necessary condition eq(16) satisfied: M full row rank => eyo E Im(M) 

-Necessary condition eq(20) satisfied 
NM_beta is full row rank (T2 does not exist) => LHS of eq(20) is 0 


[Gil G12] .ne. 0, use SQP routine in Optimization Toolbox 


>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb ,vub , "congrad" , options . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 1 0.107 (r/s) UNC = [ 0.002] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 2 0.199 (r/s) UNC = [ 0.003] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 3 0.291 (r/s) UNC = [ 0.004] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GDMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 4 0.383 (r/s) UNC = [ 0.008] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GDMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 5 0.476 (r/s) UNC = [ 0.006] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 6 0.568 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 7 0.660 (r/s) UNC = [ 0.006] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GDMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 8 0.752 (r/s) UNC = [ 0.003] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 9 0.844 (r/s) UNC = [ 0.003] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,GDMEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 10 0.936 (r/s) UNC = [ 0.005] 


>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GDMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 200 18.423 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,G0MEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 201 18.515 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
>> xio_ii ,etao_ii ,EF_ii,GDMEG_ii ,Wunc_ii ,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs ,nrcs) 
FREQ 202 18.607 (r/s) UNC = [ 0.001] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad" ,x0 ,[],[],[],[], vlb , vub , "congrad" , options ,.. . 
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>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 203 18.699 (r/s) UNC = [ 0.000] 

>> [x,FVAL,EXITFLAG] = fmincon( "fungrad 11 ,x0 ,[],[],[],[] ,vlb , vub , "congrad 11 , options ,.. . 
>> xio_ii,etao_ii,EF_ii,GOMEG_ii,Wunc_ii,bo,blkf ,ncolM,minrankNMb,minrankM,nrrs,nrcs) 
FREQ 204 18.791 (r/s) UNC = [ 0.001] 

>> [Wid_z] = ovbndunc (unc_f rd.blkf ,omeg,iv_low,iv_hig,Ts) 

Curve fit uncertainty bound # 1, blkf = [-2,0] .. 

Curve fit uncertainty bound # 2, blkf = [-2,0] .. 

Curve fit uncertainty bound # 3, blkf = [-2,0] .. 

Curve fit uncertainty bound # 4, blkf = [-2,0] .. 

Curve fit uncertainty bound # 5, blkf = [-2,0] .. 

Curve fit uncertainty bound # 6, blkf = [-2,0] .. 

Curve fit uncertainty bound # 7, blkf = [3,3] .. 

Plot summary of weights? (y/n) n 

Figure 18 shows the identified smallest additive uncertainty (circle line) subject to eigenvalue 
uncertainty allowance (of .002 on the first three structural modes) as compared to the true 
additive error (dot line). A fitted additive uncertainty model is also shown (dashdot line) 
as a reference. The level of the identified smallest additive uncertainty is very small as 
expected than the true additive uncertainty because of the eigenvalue uncertainty allowance 
of .002 on the first three structural modes was comparable to the true eigenvalue error which 
dominated the additive errors at the lower frequencies. 


Case 3 



Frequency (rad/sec) 


Figure 18: Identified smallest additive uncertainty (circle) with eigenvalue uncertainty allowance 
(.002), true additive uncertainty (dot) 
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4.3.4 Smallest additive uncertainty with uncertain mode 3 only 

As a comparison, Figure 19 shows what happens to the identified smallest additive uncer- 
tainty if the eigenvalue uncertianty allowance is increased two fold to .004 but only mode 
# 3 (resonant frequency approximately 4.5 rad/sec) as assumed uncertain with the same 
output noise allowance. Figure 19 shows that the sharp error peak near 4.5 rad/sec due to 


Case 4 



Figure 19: Identified smallest additive uncertainty (circle) with eigenvalue uncertainty allowance 
(.004) only on mode # 3, true additive Thisuncertainty (dot) 

the parametric error has been accounted for by the uncertainty allowance of .004 at mode # 
3. The other major error peak caused the parametric error near .95 rad/sec is covered by the 
additive uncertainty because apparently the large (more than twice of true parameteric er- 
ror) neighboring eigenvalue allowance did not reduce the additive uncertainty at .95 rad/sec. 
This behavior is expected since structural resonances due to distinct modes are generally 
uncoupled. Note that doubling the parametric allowance did not noticably decrease the 
additive uncertainty. Of course from a robust control design point of view, this apparently 
unnecessarily large parametric uncertainty level limits the attainable performance robust- 
ness unnecessarily. This example again demonstrates the highly non-unique nature of model 
validating uncertainty sets as highlighted in this study. 
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5 


Application Example 2: Performance validation of a ducted fan 


5.1 Caltech Ducted Fan 

In this section we outline uncertainty bounds obtained from measured data for the purpose of 
predicting the dynamic behavior of a realistic physical system, i.e. validate a system’s closed 
loop performance to within a reasonable uncertainty model for an independent controller. 
Discrepencies between responses predicted by analytical and identified models and actual 
measurements are highlighted. The predictive capability of various model validating sets 
of uncertianties are evaluated using simulated and experimental data. Despite the presence 
of a high level of aerodynamic disturbance with limited data, which significantly limited 
both identification and model validation results, the model validating uncertainties led to 
improved predictions of closed loop performance. 

A schematic of the Ducted Fan is described in Figure 20. Details on the description of the 
testbed are given in [22], Briefly, a set of encoders measure z and 9 , and a filter is used in 





z 


Figure 20: Ducted fan testbed schematic. 

the ^ channel to estimate its velocity. The measured outputs are given by ^ (airspeed), z 
(altitude), and 6 (angle of attack), about trim. A derivative filter used to estimate airspeed 
inevitably introduces a small phase lag which we can neglect in the nominal model but is 
implicitly accounted for in the uncertainty model. The inputs in the following experiments 
consist of the fan motor voltage, V m , and the paddle angle, 5 p . 
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PD Controller Inner Loop 



Figure 21: Motion about a trim at level flight stabilized by a PD controller. 


The base sampling rate for the experiments are 100 Hz which is sufficient to mitigate alias- 
ing effects with the implementation of a digital anti-aliasing filter having a break frequency 
of 5 Hz. A corresponding small phase lag was deemed a good tradeoff to mitigate allias- 
ing effects in the feedback loop during controller implementation at 20 Hz. This lower rate 
also helps to maximize the time length of recordable data, given memory limitations during 
testing while satisfying effective bandwidth requirements. Since the open loop dynamics of 
the Ducted fan is marginally stable for this particular level flight, a Proportional-Derivative 
(PD) controller is designed and implemented (for details see Appendix) in all subsequent 
experiments to stabilize and regulate its attitude, altitude, and rotation rate about a trim 
whose linearized motion is shown in Figure 21. 

Three independent sets of experiments each 20 minutes duration were completed. The 
first set involves zero test inputs and are used for identifying noise spectra in closed loop. 
The second set is used for system identification and model validation. The third set of exper- 
iments involve an additional controller over the PD controller for the purpose of performance 
prediction and validation. 


5.2 Nominal model and uncertainty structure 
5.2.1 Analytical model 

Among many choices in coordinate systems and state variables, consider an analytical model 
described by a set of first order ordinary differential equations having states, [V, z, z,9,9\, 
with inputs u T = [V r n ,5 p ] and outputs y 1 = [di/j/dt,z,9], A trim condition at level flight is 
obtained at a forward flight velocity, V a = 6 m/s, zero altitude, z, and elevator deflection, 
<5 e , and a perturbed linear model about the trim is defined. For more details, please refer to 
Appendix. 

The accuracy of an analytical model of the Ducted fan depends on many factors. A 
partial list includes the following: proper application of physical principles and selection 
of associated physical parameters governing aerodynamics on the boom, wing, and shroud, 
motor velocity; thrust dependence on applied voltage and paddle deflections, airspeeds; 
structural flexure in boom, cables, and gears; nature of friction or energy dissipation between 
moving parts including stiction; fan motor dynamics; wall effects on the flow field, etc. 
Component wind tunnel test data is used to develop aerodynamic models which is also 
not totally accurate. In addition, a difficult but crucial task is to develop a model of the 
transient aerodynamic disturbance/noise effects, which appears to account for much of the 
noisy closed loop system response. By periodically flying through its own wake in level 
flight, the resulting unsteady aerodynamic forces will depend on past control input histories, 
current attitude and position. This periodicity is evident from output noise periodograms 
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Variable 

Analytical model 

Identified model 

Vm 

0.4665 

0.4178 

8 p 

0.6745 

0.3971 

di/j/dt 

0.8341 

0.8391 

z 

0.9595 

0.6600 

e 

0.8542 

0.5156 


Table 2: Prediction errors based on I 2 norms over .1 to 10 rad/sec. 


which show peaks whose frequency matched the nominal z-axis rotation rate of the Ducted 
fan (see Figure 24). This means that the exogenous output noise is actually correlated to 
both input and output variables which will cause unknown bias errors in the identification 
of an empirical model. 


5.2.2 Identified models 

In all subsequent identification results, we apply the Subspace Model Identification [25, 15] 
algorithm. The order indicator obtained directly from time domain data show no significant 
gap to indicate a natural order of the system based on measured data in Figure 52, in 
stark contrast to the presence of a clear order for simulated data in Figure 43. Nevertheless, 
inspired by analytical models having 5 plant states plus 4 PD controller states (see Appendix 
for details), an 8th order reduced model is chosen for the system consisting of plant and the 
PD controller. As an indication of the degree of reliability in this particular identification 
algorithm, it was found that, for example, if a 9th or larger order is selected, unstable closed 
loop systems with unrealistic poles and zeros at unpredictable locations can appear. 


5.2.3 Comparison of analytical and identified models 

Figure 22 shows an identified frequency response across plant from S p to 6 with PD controller. 
The frequency response of these identified models roughly resemble corresponding analytical 
models over the bandwidth of interest, namely, .1 to 10 rad/sec. For example, at .2 and 2 
rad/sec, where the gains match well, there is about 40 to 50 degrees of phase discrepency. 
The match appears to be better in the second input, namely, paddle angle S p , than for the 
fan motor voltage input V m (see Figure 53 for the rest of the comparison). This better match 
in the second input is consistent with physical intuition, an indication of a better analytical 
model in the dynamics associated with the paddle angle input, particularly in the prediction 
of the altitude z and angle of attack 6. 

Following the notion [5] that a model is a means to predict the output response due to 
an arbitrary input, the above discrepencies do not necessarily indicate that either model 
is more accurate than the other since both predictions, of frequency responses from the 
command inputs, may differ significantly from the measured spectra. Accordingly, Table 2 
shows prediction errors of all input and output signals based on analytical and identified 
models when compared to measured signals. These prediction errors are actually I 2 norms 
over the frequency range of interest, .1 to 10 rad/sec (equivalent to a weighted RMS), and are 
normalized by the l 2 norms of the corresponding measured signals. Although the prediction 
errors for the identified model is significantly less than that of the analytical model (except 
the di^/dt channel), still there are significant discrepencies for both models when compared 
to measured signals. Since the predicted signals do not include the effects of unknown 
exogenous signals, the above discrepencies are due to an unknown mixture of model errors 
and unaccounted unknown exogenous signals, which is estimated next. 
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Figure 22: Frequency response from 5 P to 8 with PD controller: identified 8th order (solid), ana- 
lytical (dash). 

5.2.4 Equivalent output noise model 

Figure 23 illustrates the noise identification problem for a closed loop system. Figure 23a 
shows the unknown output noise q within the inner loop PD controller, K . The end goal of 
noise identification in this study, is to obtain a suitable filter, V, such that some unknown 
but unity bounded discrete-time white noise, u, sequentially drive the filters V, and plant 
output sensitivity (I — PK ) _1 , resulting in an “equivalent” output noise model suitable 
for subsequent model and performance validations and robust control application with an 
outer loop controller, C. Note that the “equivalent” output noise is intended to describe 
the particular approach whereby all unknown exogenous signals are modeled as an additive 
noise at the plant output. 

A simple (but perhaps overly simple) approach to generate an output noise model is to 
directly fit the periodograms of individual output channels under zero command input. This 
implies that the equivalent output noise model, (/ — PK)^ l V ^ is assumed to be diagonal 
which of course is not true because we know that the PD closed loop response due to 
output tracking command is strongly coupled meaning that the output sensitivity filter is 
not diagonal. An alternative and more physically compatible approach followed here is to 
assume that q is an independent set of exogenous inputs. To estimate the output spectrum 
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Figure 23: Noise model identification schematic: (a) unknown disturbance q with known inner loop 
controller K , (b) additive noise at output with inner loop closed, (c) system with an outer loop 
controller C with normalized noise v. 

of q, we utilize previously identified frequency responses of the input sensitivity and closed 
loop frequency response across plant. As outlined in the Appendix, an overdetermined set of 
loop conditions given by equations 32, are used to solve for q directly. Finally, the diagonal 
filter V is obtained by the unit norm condition of signal u. 

Figure 24 shows the identified equivalent output noise spectrum, q, based on zero inputs, 
r = 02xi, closed loop experimental data. The figure also shows the fitted diagonal normal- 
izing filter V for each channel. The identified spectra, q, indicates that there exists a large 
disturbances at about 2.5 n rad/sec where n — 1, . . . and most prominent in airspeed, dip/dt 
followed by angle of attack, 9, but not noticable for altitude, z. The frequency of these 
peaks match the rotation rate about the z-axis. Interestingly, the above noticable absence 
of periodicity in altitude is intuitively consistent. 

Finally, notice the high levels of noise (corresponding time signals are shown in Figure 
41) which is consistent with independently measured signal to noise ratios discussed in more 
detail in Appendix. The interested reader is referred to the Appendix sections for more 
discussions concerning the results and algorithms used to determine output noise models, 
analytical and identified models. The algorithms used were successfully validated on simu- 
lated data before application to actual experimental data. 


5.2.5 Uncertainty structure 

Figure 25 shows the interconnections of the uncertainty model and the equivalent output 
noise model to the inner loop PD controller, K, and the outer loop controller, C. The 
uncertainty model chosen consists of a multiplicative uncertainty at the output of the inner 
loop PD controller. A bounded but unknown equivalent noise, u, is chosen as additive at the 
output of the inner loop PD controller and is filtered by output sensitivity and normalizing 
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Frequency (rad/sec) 


Figure 24: Identified output noise q (solid), fitted noise V (dash), condition number(dashdot), 4096 
points DFT, 5 samples. 


diagonal filter V. The nominal model consists of the transfer function matrix from command 
input r to the closed inner loop output of the plant y, i.e., (/ — PK)~ 1 P. The simplified 
uncertainty structure chosen along with a zero input disturbance allowance gives the following 
augmented plant: 


V 

y 



[0,0] (I-PK)~ l P 
[0, V n ] (. I-PK)~ l P 


G(P,K,V) 



where V n := (/— P/l ) _1 V. Notice from the above augmented plant that because [Gn, G*i 2 ] = 
0, the feasibility test for model validation of this particular connection of uncertainty struc- 
ture and noise model is convex which greatly simplifies the subsequent numerical optimization 
problem. Notice that this sufficient condition is easily violated, if for example an equivalent 
input noise is defined at the input, to go with the output multiplicative uncertainty, i.e., 
Gl2 7 1 0 - 


65 






V 


J 



Figure 25: Uncertainty model interconnection about a level flight trim for Ducted fan. 

5.3 Smallest unmodeled dynamics for model validation 

Consider the problem of finding a smallest model validating output multiplicative uncertainty 
given an output noise allowance, consisting of the product of the fitted noise filter, V, shown 
in Figure 24, and an output sensitivity filter. In this study, an output sensitivity filter is 
constructed by forming the product, / 3 X 3 + T yr K where T yr is an identified model of the 
inner loop across plant and K is the given inner loop PD controller (for more details please 
see Appendix). As discussed earlier, this question can be posed as an optimization problem 
involving a linear cost function subject to a set of linear matrix inequality constraints at 
each frequency: 


min c T z 

Z 



subject to 

Qi* - IIM 2 

SiZ 

sym 

I 


\ 

sym 


Lz I \ 

where z := [Re(^); Im(^); Re(0); Im(^); x 2 ] G 7 ^ 2 ^+ 2 ^+^ anc } ^he se t Q f cons tants c, Qi, Si, 
and L are defined in equation 22, b Q is defined in equation 9, and £ 0: , is the ith uncertainty 
component of £ 0 which is defined in equation 11. For unstructured uncertainty case, r — 1, 
while for diagonal uncertainty case, r = 3. 

The interior-point algorithm found in MATLAB’s LMI Toolbox [20] is used to solve for 
the global minimum. With multiplicative uncertainties for all output channels, the existence 
of a model validating set is obviously not in question and the interest lies in determining the 
necessary magnitude level and frequency shape of the unmodeled dynamics (which implicitly 
includes some unknown level of contribution from parametric errors in the nominal model). 
Table 3 show the various cases considered which arises from a combination of nominal model, 
assumed uncertainty structure, and a choice on the level of output noise allowance based on 
a factor of fitted noise spectrum. 
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Case 

Nominal 

Uncertainty str 

Noise 

Expl-2-D 

Expl-3-D 

Expl-4-D 

Analytical 

Analytical 

Analytical 

Diagonal 

Diagonal 

Diagonal 

,5x fitted 
fitted 
2x fitted 

Exp4-2-D 

Exp4-3-D 

Exp4-4-D 

Identified 

Identified 

Identified 

Diagonal 

Diagonal 

Diagonal 

,5x fitted 
fitted 
2x fitted 

Expl-2-F 

Expl-3-F 

Expl-4-F 

Analytical 

Analytical 

Analytical 

Unstructured 

Unstructured 

Unstructured 

,5x fitted 
fitted 
2x fitted 

Exp4-2-F 

Exp4-3-F 

Exp4-4-F 

Identified 

Identified 

Identified 

Unstructured 

Unstructured 

Unstructured 

,5x fitted 
fitted 
2x fitted 


Table 3: Cases considered for computing smallest output multiplicative uncertainty. 

Figure 26 shows the smallest model validating diagonal output multiplicative uncertainties 
based on analytical and identihed nominal models. Assuming unstructured output multiplica- 
tive uncertainties instead, Figure 27 shows the smallest model validating maximum singular 
values of the uncertainties based on analytical and identihed nominal models. The effects 
of using different levels of equivalent output noise allowances are also shown: the diamond 
solid line denotes uncertainties obtained based on a unit factor of fitted noise allowance, the 
triangle dotted line is based on twice the fitted noise allowance, and the square dotted line 
is based on half the fitted noise allowance (refer Table 3). Based on the results in Figures 
26 and 27 we make the following observations: 

1. Larger noise allowance leads to smaller model validating unmodeled dynamics. How- 
ever, at twice the fitted noise allowance, the identihed nominal model case gives unreal- 
istically small levels of unmodeled dynamics (say 1 % or less at many frequency points). 
This also means that the validation was likely achieved mostly with noise allowance, 
which could lead to unrealistic predictions of robust stability. 

2. Identihed nominal model cases resulted in smaller uncertainty levels than analytical 
nominal model, for both types of uncertainty structure and three different noise al- 
lowance levels. Since the uncertainty ball of the identihed nominal model is smaller 
than that of analytical nominal model, the identihed nominal model appears to be more 
accurate than the analytical nominal model. This result is consistent with the earlier 
prediction error analysis of nominal models. 

3. The maximum singular value of the unstructured uncertainty cases resulted in smaller 
uncertainty levels than structured uncertainty cases (its maximum over all channels). 
This is not surprising since unstructured uncertainty has a larger number of uncertainty 
components than diagonal uncertainty. However, a smaller normed uncertainty set do 
not necessarily imply a smaller set of plants if the number of uncertainty components 
are different. 

4. At lower frequencies, say < 2 rad/sec, smaller uncertainty levels are required for both 
types of nominal model and uncertainty structure. This may be due to the following 
reasons: 

• The frequency responses of both analytical and identihed nominal models have 
larger gains at lower frequencies (see Figures 53), and the multiplicative uncer- 
tainty at output represents a factor of these responses. 
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• The nominal model is actually more accurate at lower frequencies; for example, 
the analytical and identified transfer functions from 5 p to 6 match better at lower 
frequencies between .1 to 2 rad/sec then over 2 to 30 rad/sec range (see Figure 
22 ). 

• The use of a larger noise allowance at lower frequencies which inadvertently can re- 
duce model uncertianty level. A larger identified noise spectra at lower frequencies 
is indicated in Figure 24. 


Exp1-2D( — ), Exp1-3D(-), Exp1-4D(:) 12000 pts Exp4-2D(— ), Exp4-3D(-), Exp4-4D(:) 12000 pts 



Figure 26: Smallest diagonal output multiplicative uncertainty about analytical (left) and identified 
(right) nominal models: noise allowances of .5x fitted (dash), fitted (solid), 2x fitted (dot). 

In summary, the above model validating uncertainty levels obtained from experimental 
data for the Ducted fan is not small. Simulation results (see Appendix for details) also 
indicate that these results can be sensitive and therefore unreliable in the presence of a high 
level of noise in data, exactly, what was observed in dealing with measurement data. This 
large uncertainty level appears to be due to two reasons: 

1. significantly inaccurate nominal models which is evident from early prediction error 
analysis of both models, and 

2. large exogenous disturbance and noise which is evident from earlier construction of 
noise models which causes DFT errors from nonperiodicity. 

We have also shown that this lack of periodicity cannot be eliminated by innovative tapered 
time windowing. In fact, simulation studies indicate that with a high level of noise in data, 
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Exp4-2F( — ), Exp4-3F(-), Exp4-4F(:) 12000 pts 



Figure 27: Smallest unstructured maximum singular value of the output multiplicative uncertainty 
about analytical (left) and identified (right) nominal models: noise allowances of .5x fitted (dash), 
fitted (solid), 2x fitted (dot). 


it is difficult to recover the correct model error, even with an accurate nominal model and 
correct uncertainty structure. On the other hand, simulation studies also suggests that 
incorrect or indeterminate uncertainty structure do not preclude satisfying model validation 
conditions. For example, assuming unstructured uncertainty when the correct model error 
is diagonal resulted in a smaller nornred unstructured uncertainty (when compared to the 
norm of true diagonal model error), exactly what was observed in dealing with measurement 
data. Nevertheless, since the important property of an uncertainty model is in its predictive 
capability of closed loop performance rather than its particular uncertainty structure or 
size, we next show the results of its predictive capability of closed loop performance for an 
independent outer loop controller. 
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5.4 Performance validation analysis 

In this section, we show comparisons between measured worst case performance and various 
predicted performance bounds based on corresponding model validating uncertainty sets. 
The measured worst case performance is defined here as the maximum possible closed loop 
amplification over bounded signals. In contrast, the predicted performance bounds are based 
on the predicted worst case signal amplification over a corresponding model validating set 
of plants. Closed loop refers to a second controller loop with C around the inner controller 
loop involving K . In the Ducted fan, the inner loop controller was necessary since the open 
loop was marginally unstable to collect data for system identification. 

Figure 28 illustrates how closed loop performance can be validated. The idea is that the 
measured worst case disturbance amplification of the closed loop system is bounded by (or 
comparable to) the predicted worst case disturbance amplification based on an uncertainty 
model. This predictive capability is an important necessary characteristic of a useful uncer- 
tainty model in the context of robust control. The predicted worst case performance for a 



Figure 28: Performance validation schematic for Ducted fan: predicted versus measured 
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given nominal and uncertainty model is illustrated in Figure 29 where skewed-/i [1] is defined 
at each frequency and is given by 

/j, s := {min p : det (I - R(p)M(P , K, V, C)A B W) = 0, d(A s ) < l} -1 (25) 

= max a[W perf P u (M(P, K, V,C), A unc W unc )] (26) 

cr(A uri c)<l 

Aunc^'D 


where 


R{p) 

M(P, K, V, C) 


0 

Plperf 


Pi 


G(P,K, V)[I,I], C 


G(P , K, V ) denotes the augmented nominal system consisting of the inner loop whereas 
M (P, K , V, C) denotes the nominal closed loop system with the second controller. Equation 
(26) indicates that p s denotes the worst case maximum singular value frequency response at 
a given frequency. This worst case is over the uncertainty model set, T>, while the maximum 
singular value reflects the multivariable nature of input and output signals and its spatial 
dependence. Equation (25) connects p s to the more well known /! which is typically used to 
compute the former by a scaling procedure described in more detail for example in [1], 



Figure 29: Computation of worst case performance from r to y (top), canonical form (middle), 
performance channel scaled and closed (bottom). 


From equation 26 and middle figure in Figure 29, note that p s includes a directional worst 
case with respect to unknown exogenous inputs, /3. For performance validation purposes, 
this is a bit difficult to experimentally measure the worst case responses with respect to 
unknown exogenous inputs. Hence, although mathematical models for the exogenous noise 
are available, we consider only the measurable worst case response from command input, r, 
to measured output, y, corresponding to the p s with performance from r to y only. Therefore, 
the following performance validation results do not include the responses due to unknown 
exogenous signals, which is known to be significant. Furthermore, since it is a non trivial 
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matter to test all possible command inputs in the hope of directly measuring the worst case 
signal amplification, we circumvent this by identifying an empirical model of the closed loop 
system with the goal of computing the maximum singular value of its empirical frequency 
response matrix, resulting in an indirect measurement of worst case directional response. 

Figure 30 shows a comparison of the measured worst case response (line), to the predicted 
responses: nominal (dash), fi s based on structured uncertainty (line-circle) and unstructured 
uncertainty (line-square). These predicted responses are based on an analytical nominal 
model. Figure 31 shows similar comparisons but its predicted responses are based on an 
identified nominal model. Both predicted responses are based on an assumed noise allowance 
of fitted noise with factor 1. Based on the results in Figures 30 and 31, we make the following 


MEASURED maxSV (-); PREDICTED: Nominal (— ), bnds struct unstr (:); Y1 X3 F ex 



Figure 30: Worst case response from r to y with noise allowance fit factor 1: analytical nominal 
model with diagonal and unstructured uncertainties. 

observations: 

1. Without using uncertainty models, the predicted performance based on identified nom- 
inal model (dash line in Figure 31) gives a much better match with measured response 
(solid line) than performance based on analytical nominal model (dash line in Fig- 
ure 30). This is consistent with both pervious analysis of uncertainty size levels and 
prediction errors. Most prominently, the analytical nominal model predicts an overly 
pessimistic response at lower frequencies and overly optimistic response at higher fre- 
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MEASURED maxSV (-); PREDICTED: Nominal (— ), n bnds struct unstr (:); Y4 X3 F ex 



Figure 31: Worst case response from r to y with noise allowance fit factor 1: identified nominal 
model with diagonal and unstructured uncertainties. 


quencies beyond 2 rad/sec. On the other hand, the identified nominal model predicts a 
slightly optimistic response at most frequencies except over a narrow bandwidth within 
1 and 2 rad/sec (see Figure 31). 

2. Using uncertainty models, the predicted worst case response (i) worsens the pessimism 
at lower frequencies but improves the overly optimistic prediction at higher frequencies 
for the analytical nominal model case (dot and dash-dot lines in Figure 30), and (ii) 
improves the slightly optimistic prediction at lower frequencies but worsens the slightly 
optimistic prediction at higher frequencies for the identified nominal model case (dot 
and dash-dot lines in Figure 31). Since errors at lower frequencies (say < 2 rad/sec) 
dominate and therefore its prediction more significant, we infer that using uncertainty 
models improved the identified nominal mode based prediction while it increased the 
conservatism of the analytical nominal model based prediction. 

3. The predicted worst case performance locally peaks and then drops significantly to 
the predicted nominal levels between 1 to 3 rad/sec. The local drop in worst case 
performance corresponds to the large resonant noise allowance (Figure 24) leading to 
smaller model validating unmodeled dynamics. Conversely, the local peaks of the worst 
case performance (Figures 30 and 31) correspond to a local drop in the noise allowance 
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(Figure 24) leading to larger model validating unmodeled dynamics. This is a result 
of treating noise as a given allowance in computing the minimum unmodeled dynamics 
necessary to satisfy model validation conditions. 

4. The deviation of predicted responses based on identified nominal model (dot and dash- 
dot lines in Figure 30) is significantly less from its corresponding nominal predicted 
response than the deviation of the predicted responses based on analyical nominal 
model (dot and dash-dot lines in Figure 31) to its corresponding nominal predicted 
response. This smaller deviations of predicted responses is consistent with the smaller 
sizes of the model validating uncertainties for the identified nominal model case. 

5. The deviation of the predicted responses based on unstructured uncertainty (dotted 
lines) from the predicted nominal response is significantly smaller than the correspond- 
ing deviation based on diagonal uncertainty (dash-dot lines) for both types of nominal 
model. This deviation is particularly large in the analytical nominal model based pre- 
dictions at lower frequencies; perhaps the assumed diagonal uncertainty structure could 
not easily account for a large coupled dynamics error that occurs for a analytical nomi- 
nal model. Based on comparisons to measured responses (solid lines), the unstructured 
uncertainty (dotted lines) appears to be more suitable for use with analyical nominal 
model while the diagonal uncertainty (dash-dot lines) appears to be more suitable for 
use with identified nominal model. 

Figures 32 to 35 show the effects of varying levels of noise allowances (.5, 1, 2) of fitted 
noise, on the closed loop performance, for both types of nominal models and both types of 
uncertainty structure. The measured worst case is shown as a reference by a solid line. 
From the four figures, note the following observations: 

1. With increasing noise allowance (from dot to dashdot to dash lines), the level of un- 
modeled dynamics decreases and the predicted worst case response decreases for both 
types of nominal models and uncertainty structures. 

2. Varying the noise allowance levels did not significantly improve the overly pessimistic 
predictions based on analytical nominal model (Figures 32 and 33) mainly because 
the level of pessimistic predictions dominates. However, note that with the use of 
smaller noise allowance (.5), the prediction based on analytical nominal model at higher 
frequencies (say > 2.5 rad/sec) matched the measurement more closely, although less 
significant. 

3. Using a noise allowance level of 1.0, the predicted response, based on identified nominal 
model with diagonal uncertainty (Figures 34), closely matched the measured response. 
In contrast, this “correctly” scaled noise allowance level of 1.0, did not appear to give 
better results if unstructured uncertainty is used (Figure 35). 

4. The predicted worst case responses based on unstructured uncertainties (Figures 33 
and 35) appears to be less sensitive to noise allowance levels for either type of nominal 
model than diagonal uncertainties (Figures 32 and 34). 

As a summary based on the above observations, the predicted worst case closed loop perfor- 
mance based on a model validating diagonal uncertainty about the identified nominal model 
with the more realistically scaled noise allowance (i.e. fit factor of 1.0), gave the closest 
match to the maximum singular value of the identified closed loop system at lower frequen- 
cies where the outputs were largest. The presence of a high level of exogenous disturbances 
due to unknown unsteady aerodynamics, concentrated near the rotation frequency of about 
2.5 rad/sec, is a significant obstacle in the uncertainty level determination methodology and 
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MEASURED: maxSV (-), PREDICTED: (i. bnds with noise levels .5(:), 1(-.), 2(— ); Y1 X2-4 D ex 



Figure 32: Worst case response from r to y for analytical nominal model with diagonal uncertainty: 
effects of noise allowances .5 (dot), 1 (dashdot), 2 (dash), and measured (solid). 


subsequent performance validation near this frequency. On the other hand, unknown dis- 
turbances always present in physical systems is always a limiting factor in model validation. 

Despite some success in the above performance predictions, the selection of a set of uncer- 
tain parameters and structure of unmodeled dynamics is not completely clear. Furthermore, 
any particular uncertainty parameter and structure that leads to a model validating set is 
likely to be redundant and there is no clear “best” selection suggested by the physics of 
the Ducted fan testbed. Nevertheless, experimental results appended by simulation results 
suggest that it may not be necessary to know the “correct” uncertainty structure, to obtain 
uncertainty models that are useful for robust control applications. 

Finally, the limitations of the results and scope of this study is noted that further research 
could address. This includes (i) performance validation is limited to worst case response 
comparisons based on maximum singular values and could be extended to gain and phase 
variations of each input-output channel, (ii) a procedure for the performance validation 
involving the worst case unknown exogenous disturbance and noise is not clear and therefore 
not considered, (iii) an improved procedure for obtaining an equivalent output noise model 
from closed loop measurements is desirable, and (iv) techniques for distinguishing exogenous 
random noise from model error effects from limited closed loop response data is desirable. 
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A Ducted Fan in level flight 


A schematic of the Ducted Fan is described in Figure 36 and Figure 37. (0, A", Y, Z ) is an iner- 
tial coordinate system, (O s , X s , Y s , Z s ) a stand center of mass, body fixed coordinate system, 
(Ob, X 'b, Yb, Zb) a shroud center of mass, body fixed coordinate system, and (O w , X w , Y w , Z w ) 
the wind coordinate system. -0 is the angle between [OX) and the projection of [O s X s ) on 



Figure 36: Ducted Fan Body 

the plane (XOY), 0, the pitch axis angle i.e. the angle between the local horizontal [O s X s ) 
and the [0&X&) axis and z the algebraic distance 00 s . The airspeed at the center of mass 

of the ducted-fan (shroud+wing+boom) is denoted by V s = — , the angle of attack by a s , 
the flight path angle by %, the paddle angle by 5 P , the elevator angle by 5 e and the motor 
voltage by V m . Subscripts s and w denote variables related to shroud and wing respectively. 
The Lagrange’s equations of motion are given by 

lzfi> = -rsF s Xs -r w Fl-r b F b Xs -dI p x n cos 9 (27) 

777 / 

(mb + nrif + — )z = mg + F s Zs + F Zs + F b Zg (28) 

4 br 9 = My s + My s + I Xb ^O. cos 6 (29) 

where I s z is the moment of inertia of the fan, boom and counterweight about OZ ; /£ , 
the moment of inertia of the fan about OfeFj, when wing in rearward position; I x , the 
moment of inertia of the propeller and the motor about OfcA^; F x , the X s component of the 
boom aerodynamic force; F Xs and F z , the X s and Z s components of the wing aerodynamic 
force; My , the Y s component of the wing aerodynamic moment; F x and F z , the X s and 
Z s components of the shroud aerodynamic force; Mp s , the Y s component of the shroud 
aerodynamic moment; r^, the effective moment arm for the boom; r w , the effective moment 
arm for the wing; r s , the distance between Ob and the plane XOZ ; mj the mass of the boom; 
m c , the mass of the counterweight; m/, the mass of the fan; r, the pulley gear ratio; g, the 
gravity and the motor velocity; m := mb Fmf — 
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Figure 37: Shroud and Wing Conventions 


Boom and wing aerodynamics Due to the rotation of the fan, the wing will experience 
increasing velocities from its root. A similar situation occurs for the boom. The aerodynamic 
forces of the wing and boom may be written as: 

F b Xs = F b Xs (iP,z,C Db ) 

FI = F$'(i,,z,C%,CZ) 

FI = Fg'(i,,z,C%,C%) 

Ml = Ml(iP,z,CZ,C™,C™) 

where Co b is the drag coefficient of the boom (cylinder) and j, C™ and C^, determined 
by a series of wind tunnel tests, are functions of a w and 6 e . For now, the effects of all other 
derivatives are assumed negligible. 


Shroud Forces and Moments There are two contributions to the shroud aerodynamic 
forces and moments: the thrust from the ducted fan engine and the aerodynamic forces of 
the shroud. In order to determine the shroud forces FL and Ff , and moment My, four 
2-D table lookup (F|T, C S D , C S L and C S M ) as a function of a s , V m or 5 P are necessary. 


Motor Speed A 2-D table lookup, function of V m and 5 P , will be used to determine the 
motor velocity D. The effects of V and a are neglected. 

The actuators are assumed to be perfect and the moments of inertias of the fan (shroud 
and wing) and the propeller about ( ObX b ) and ( ObZ b ) are neglected. 
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A.l Linearized, unstable model about trim 


A ducted fan analytical model can be described from equations 27 to 29 by a set of first 
order ODEs x = f(x,u ) where x T = [V, z, z, 9, 6] and u T = \V mi 5 P , 5 e ] where V := ipr s 
and r s = 2.35 is the radius in meters to convert angular to linear velocity. Although the 
horizontal linear distance navigation variable, ipr s which is the integral of V, is not necessary 
as a state, the navigation variable, z, is necessary to control the altitude during experiments 
since the vertical movement is constrained by the testbed. Figure 21 is a schematic of the 
linearization about a trim point. A trim condition at level flight is obtained by solving 
the above non-linear algebraic equations at a forward flight velocity, V a = 6 m/s, and zero 
altitude, z, and elevator deflection, 8 e , for the pitch angle 9, the motor voltage V m , and the 
paddle angle S p . In other words, the trim conditions satisfy the following conditions: 


trim [Vo, 0; 0) ^trimi 0] > 


U 


trim 


= [W. 

rim "> ^ Ptrim "> ^ e trim ] , f Ot rimi ^ trim ) 0lx5 


By computing the Jacobian of / about the trim point a perturbed linear model was 
obtained about, V = 6 m/s, z — 0 m, 6 — 9.8104 deg, V m = 1.0927 Volt, 8 P = —9.1725 deg, 
S e = 0 deg, using (linmod.m) with the following results: 
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V M 


0.5569 -0.0088 

0 0 

-0.0344 -0.1432 
0 0 

0.7034 -2.3021 



Table 4 show the eigenvalues and eigenvectors of the linear continuous time plant model. 
This model predicts the following: 

Ay rigid body mode in z-axis 

A 2 , 3 : (stable, unstable) pair of slow (uj — .02 Hz) z-ip axes coupled mode; Phugoid motion 

Az^: damped oscillating (a; = .15 Hz,( = .35) 9-z-ip axes coupled mode 



Ai 

A 2 

A 3 

A4,5 

Degree-of-Freedom 

0 

-0.1195 

0.1370 

-0.3507 ± 0.9127i 

5V 

0 

0.0989 

0.1060 

-0.0637 - 0.11771 

Sz 

1 

0.9880 

0.9846 

-0.6390 + 0.2532i 

Sz 

0 

-0.1181 

0.1349 

-0.0070 - 0.67211 

se 

0 

0.010 

-0.0332 

-0.1625 + 0.056U 

se 

0 

-0.0013 

-0.0045 

0.0058 - 0.1680i 


Table 4: Eigenvalues and eigenvectors of the unstable, continuous, analytical model. 


The encoders directly measure "0, z and 6, and a filter is used in the ^ channel to estimate 
its velocity. Hence the measured outputs are assumed to be z and 9, about trim. This 
derivative filter inevitably introduces a phase lag which we ignore for simplicity. In addition, 
the encoder outputs and actuator inputs are sampled and held at 100 Hz. This sampling rate 
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appears sufficiently fast to mitigate aliasing effects with the implementation of a digital anti- 
aliasing filter having a break frequency of 5 ffz. The filter was designed after preliminary 
testing and analysis of an analytical model of the Ducted fan which indicated that the 
system’s significant dynamics was limited to about 2 ffz and it was decided that a discrete 
time controller will be designed and implemented at 20 ffz. A small phase lag introduced 
by this filter was deemed a good tradeoff to mitigate alliasing effects in the feedback loop 
during controller implementation. 

Finally, it turns out that the above trim point is marginally unstable (or stable ?) and a 
stabilizing controller is required to collect any reasonable length of measurement. 


A. 2 Truth plant model and controller for simulation 


The previous unstable analytical plant model is taken as the unknown true plant for gen- 
erating simulated data. A truth plant model is defined for the purpose of simulation study 
only and is intended to define a particular idealized physical system and is assumed to be 
uniquely defined by a given mathematical model. 1 The main goaf in this section is to deter- 
mine how well (and validate) system identification algorithms recover a given true plant from 
simulated measurement data, before we consider actual measurement data. Although the 
physical states in general cannot be recovered through system identification, we identify and 
compare coordinate invariant physical properties such as plant eigenvalues and input-output 
maps such as frequency response. In addition, since the measured outputs are physical 
variables, a physical interpretation of the mode shape is possible. 

The truth plant model is discretized at 20 Hz, which corresponds to the sampling rate 
of an actual laboratory experiment. The following Proportional-Derivative (PD) controller 
discretized at 20 Hz is used in the generation of the subsequent closed loop simulation data. 
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This controller consists of two separate feedback loops, namely, the first input SV m feeds 
back only Sip while the second input 5 p feeds back a linear combination of Sz and 59. The first 
input is designed to regulate the constant rate rotation by providing a negative proportional 
feedback on the deviation of the flight velocity about trim while the second input is designed 
to stabilize and regulate both its attitude and altitude. It is physically significant to note 
that this particular controller actually stabilizes the Ducted fan in a laboratory experiment. 
The closed loop performance of this PD controller is used as a baseline to validate the 
performance of an independent controller in an outer loop. 

4 However, in many physical systems which are of engineering interest, a “true” plant defined by a unique math- 
ematical model may not exist. Hence, much debate revolves around the mismatch between measured data and that 
predicted by various “true” mathematical models, which motivates voulmnous previous work in areas such as system 
identification, model validation, robust control and countless related experiments. 
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A. 3 Comparison of analytical model response to measurement 

To get an idea of the accuracy of the analytical model, we compare the spectra of the 
measured closed loop output (input) to the predicted closed loop output (input) from the 
analytical model driven by the same test signal used in the experiment. The unstable 
analytical plant model with the PD controller used in the experiment is stable in closed 
loop. No unknown noise or disturbance is assumed in the simulated response. 

Figures 38 and 39 show measured (solid) and predicted (dash) magnitude spectra of the 
closed loop outputs and inputs, respectively. Specifically, note that 


measured(solid), analytical(dash), difference(dot), 4096 points @ 20Hz 




Figure 38: Measured magnitude spectra of PD closed loop output vs analytical model prediction 


• The first channel in the measured output, 5dif}/dt, and the fan voltage input, V m , show 
resonant peaks at about 2.5 rad/sec which corresponds to the nominal rotation rate of 
the Ducted fan about the z-axis, i.e., difj/dt = V/r s = (6m/s)/(2.35m) = 2.55 rad/sec. 
There also appears to be higher harmonic resonant peaks at about 5 and 10 rad/sec. 

• The discrepencies, as a percentage of its nominal response in the output spectra, is 
generally larger with increasing frequencies. This is likely due to the combination of (i) 
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Figure 39: Measured magnitude spectra of PD closed loop inputs vs analytical model predictiction 


a less reliable analytical model at higher frequencies and (ii) the presence of wideband 
noise which plays a more prominent role at higher frequencies where the plant rolls off. 


In summary, there are significant discrepencies between the spectra of the measured closed 
loop output and the output simulated from an analytical model in all five channels. Note 
that similarities in the spectra magnitudes may appear more optimistic than it actually is 
since their phases could be totally off. Hence, the figures also show the magnitude of the 
differences in their complex spectra (dot). To get some feel for model errors from these 
signals, we compute the following signal ratios 


| Vmeas Vs 
1 1 2/szjti 1 1 
| V mens Vs 


mult ) 
^"(^add) 


which are given in Figures 40. The above ratios roughly approximates a lower bound on the 
2-norm of a complex 3 by 3 full-block output multiplicative uncertainty and the 2-norm of a 3 


84 


levl/lvl (solid), levl/lrl (dash), max SV Diant (dot), 4096 ooints & 20Hz 



Frequency (rad/sec) 


Figure 40: Ouput difference signal ratios. 


by 2 additive uncertainty, respectively. The figure shows that the multiplicative uncertainty 
increases with frequency while in the additive uncertainty it decreases. This makes sense 
because the multiplicative uncertainty is a factor of the noiminal plant which rolls off as 
indicated by the maximum singular value frequency response of the closed loop system using 
analytical plant model (dotted line). 

To get a feel for the noise and disturbance levels in the system, Figure 41 shows the 
measured closed loop time histories when the system is forced (left column) and free (right 
column). Clearly, the system is very noisy particularly in the forward motion (first channel). 
The perturbed forward motion (dd^/dt), in both the forced and free closed loops matche the 
rotation rate about the z-axis. This is more evident from an estimated equivalent output 
noise spectrum which show peaks whose frequency matches the z-axis rotation rate of the 
Ducted fan. By flying through its own wake periodically, the resulting turbulent aerodynamic 
forces will depend on past control input histories (fan speed and vector thrust angle) and 
current attitude and position. This causes the exogenous output noise to be correlated to 
both input and output variables which will cause unknown bias errors in the identification 
of an empirical model. In addition, a vertical drift or a very slow oscillation is also present 
in the z plot of free response. 

In summary, a large discrepency exists between the predicted spectra from the analytical 
model and measured data which suggests a significant level of inaccuracy in the analytical 
model. Furthermore, the unknown noise and disturbances in the system appears significant 
indicating that an empirical nominal model development through system identification may 
not be trivial, and properly accounting for unknown exogenous disturbances in any validation 
process is clearly important. Of course, the accuracy of any analytical model of the Ducted 
fan depends on the assumptions made about the physical system during model development. 
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Figure 41: Forced (left column) and free (right column) closed loop output measurements. 

A partial list of these assumptions involve: 

• inertial properties, dimensions, and center of gravity location, . . . 

• nature of friction, energy dissipation mechanism, stiction, between rotating parts, . . . 

• structural flexure in the composite boom, suporting cables, gears, . . . 

• aerodynamic models 

— boom aerodynamics, F\ ( / ip,z,CD b ) 

~ wing aerodynamics, F^, z,C%,Ct), (^, i, Q) 

— shroud aerodynamics, Fy , F| , My 

• motor velocity dependence, £l(V m , S p ) 

• wall and turbulent wake effects on aerodynamics 

• uncertainties in tabulated component wind tunnel test data 

Therefore, the Ducted fan is difficult to model accurately with confidence because it is 
difficult or near impossible to verify all the assumptions related to the above list. 
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A. 4 Identified models from simulated data 


In this section, we investigate whether an empirical model obtained directly from experimen- 
tal data can reduce the significant discrepency that exists for the analytical model. We show 
system identification results based on simulated data to validate the identification method 
which is subsequently applied to experimental data. Based on both simulation and exper- 
imental results, we discuss inherent limitations in the system identification approach. All 
closed loop simulated data in this chapter is generated from an assumed unstable analytical 
truth model with a known PD control law. 


A.4.1 System identification algorithm 

We briefly outline a popular a time domain identification method called the Subspace Model 
Identification (SMI) [25, 26] which we use to obtain all subsequent system identification 
results in this report. The identification software SMI Toolbox Version 1.0 was obtained 
from [27]. The SMI method is based on approximating structured subspaces of observability 
matrix from input-output data and then directly realizing state space matrices from the data 
equation with the aid of a shift-invariant property of the observability matrix. An output- 
error state space structure is assumed as shown in Figure 42 which also covers black-box 
models such as ARMAX and OE [5] . The following is a description of the ordinary MOESP 


Input 


zero-mean white noise 


Noise Shaping 
Filter 


colored noise 


Unknown 

Deterministic 

System 




Output 


Figure 42: Output Error State Space model identification (MOESP) setup, 
algorithm in a nutshell: 

1. Suppose we have measured input and output data sequence from time index j to 
(j + N + i - 2) 

[Uji • • • j 'U j j+N-\-i—2\ ? \]Jj 5 • • • j Vj+N+i— 2 ] U'k ^ R Vk ^ R 

2. Assume that this finite segment of the data sequence satisfy the noise- free LTI relation 

Assumption 1: x^+i = Ax fc + Bu y k = Cxk + Duk 
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meaning that the following data equation holds: 
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(30) 


3. To obtain an approximation to column space of T*, assume the following: 


Assumption 2: 
Assumption 3: 

Assumption 4: 


(A, B,C, D) is minimal of order n 
i > n 


rank 


U Ui , N 
A A,n 


mi + n 


and a RQ factorization of the fat (by assuming a sufficiently long measured time seg- 
ment, i.e., a large N ) data matrix pair is computed 


Ui^N 


i?H 
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Q 1 

bl,j,AT 


R 2 I 

R 22 


Q 2 


where Rn and R 22 are lower triangular matrices and the rows of Q 1 and Q 2 are or- 
thonormal. With Assumptions 1 to 4, it can be shown that Image(R 22 ) = Image{Ti ) 
and an orthonormal bases for Image{Ti ) can be computed from the SVD 


R 22 = [ U n U ± } 
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4. 


The orthonormal base matrix, U n , which is related to the unknown observability matrix, 
Tj by U n = TjT where T is an unknown square nonsingular matrix, can be used to com- 
pute state and output matrices A T and Ct using a structural property of observability 
matrix 
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(31) 
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Since T is square and nonsingular, the above “shift-invariant” property can be rewritten 



so that At can be computed uniquely since Un has full column rank per assumptions 
2 and 3. The corresponding output matrix, Ct , can be determined by noting that U n 
is the observability matrix in transformed coordinates and therefore the partitioned 
matrix consisting of the first /-rows of U n is Ct- 

5. Notice from the noise-free data equation (30) that if A and C are known, the only 
unknowns are B and D which appears linearly. Hence, after some algebraic gymnastics 
[25], one can obtain the following overdetermined equation ((Zi — n)i x m equations 
with {m + n) x m unknowns) 
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where 5 := {U^) and a least squares solution of B T and D can be obtained. 

It is important to note that a least squares solution will result in satisfying only ap- 
proximately both the shift invariant property, equation (31), and the noise-free data 
equation (30) and this limitations on model accuracy is unclear. 

Overall, SMI approach is computationally easy but the performance of the method is un- 
proven under noise and disturbance (SMI is a deterministic approach), nonperiodic test 
signal, uncertainty in the assumed minimal order of the state. Various recent work reported 
in IFAC’s Symposium on System Identification SYSID 2000 [15] is a testament to the high 
level of interest in this particular identification method. Hence, before we consider experi- 
mental data, we evaluate this algorithm using simulated data. 


A. 4. 2 Identified models of (/ — KP) 1 from simulated data 

Simulated data containing low (variance cr = 10~ 8 ) and high (variance a = .03) measurement 
noise levels are considered. The simulated high noise level roughly approximates the actual 
measurement noise levels in terms of signal to noise ratio. In the system identification 
simulations for the high noise level cases, 100 minutes of data (or 29 signal samples each 
approximately 3.4 minutes duration sampled at 20 Hertz) are assumed while for the low noise 
level cases only 20 minutes of data (or 5 signal samples) are assumed. In actual experiments, 
the measurement data is 20 minutes duration. 

Figure 43 shows the singular values which help determine suitable orders for the input 
sensitivity state space model corresponding to low (left figure) and high (right) noise cases. 
In both cases, a (carefully selected at) random upper bound order of 18 is chosen to compute 
the singular values of the Hankel matrix. Notice from the low noise case (left figure) that 
nine singular values corresponding to the true closed loop system states stands out because 
of the very low noise floor. This is in sharp contrast to the high noise case (right figure) 
where the high noise floor conceal states beyond the first five. 

In the low noise case, it was found that although the largest singular value gap occurs after 
the first nine states, thus suggesting a natural reduced order, both 9 and 8th order realizations 
were unstable and the largest stable realization order was 7. In the high noise case however, 
the largest singular value gap occurs after the first five states but this gap is relatively small. 
Therefore in both cases, all 9 states of the truth model was not recovered. In retrospect, 
this limitation is expected in most nontrivial applications in system identification where the 
lesser observable and/or controllable states are difficult to recover due to their relative lack 
of participation in the output signal. 

Figures 44 and 45 show the frequency responses for the true and identified 7th order 
input sensitivity model corresponding to low and high noise cases respectively. The dash 
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System order 



System order 


Figure 43: Singular values to determine order of input sensitivity from simulated data, low noise 
a = 10” 8 (left), high noise a = .03 (right). 


lines denote the true frequency response while the solid lines denote the identified frequency 
responses. It is clear from these figures that identification results based on low noise data 
give a much more accurate model. The gains of the difference between true and identified 
frequency response for the low noise case is about one percent of the true frequency response 
(analogous to one percent multiplicative error for each channel) over all frequency. On the 
other hand, for the high noise case, the match at the low frequencies are poor while the 
match at the higher frequencies for the diagonal channels are excellent. These appear to 
be the result of the input sensitivity which attenuates the lower frequency signals at the 
input, causing a poor signal-to-noise ratio and the characteristic loop gain rolloff which 
drives the input sensitivity to identity, causing a good signal-to-noise ratio. Clearly, inspite 
of a significantly longer data record, using data with a high level of noise (comparable to 
measured data), resulted in a significantly less accurate model as compared to resulting 
model based on data with a low level of noise. 


A. 4. 3 Identified models of (/ — PK ) 1 P from simulated data 

Figure 46 shows the singular values which help determine the orders of a state space model of 
the closed loop transfer function across plant based on simulated data containing low (left) 
and high (right) noise levels. Again, an upper bound order of 18 is chosen to compute these 
singular values. For the low noise case (left figure) 10 states stands out (true states being 9 
only) while only 5 states stands out . for the high noise case where the noise floor is high. For 
the low noise case, it was found that 7, 8 and 9th order realizations were unstable and so a 
6th order stable realization was chosen. On the other hand, a stable 7th order realization was 
chosen for the high noise case although it appears that only five states appear to dominate 
the response. In general, choosing higher order models often introduced fictitious poles and 
zeros and did not improve the model. 

Figures 47 and 48 show the true (dash) and identified (solid) frequency responses for the 
6th and 7th order models of the closed loop transfer function across plant corresponding to 
low and high noise cases respectively. Similar to earlier results for input sensitivity, low noise 
data leads to a more accurate model (inspite of using much longer data record for high noise 
data case). The multiplicative error is roughly 1 percent error for the low noise case whereas 
it is roughly 10 % error for the high noise case, in the frequency range of interest (say .1 to 
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Figure 44: Input sensitivity, 20 minutes simulated data: True (dash), Identified 7th order (solid), 
low noise cr = 10~ 8 . 

10 rad/sec). 


A. 4. 4 Summary of simulation results 

The results based on simulated data appears to validate the SMI identification methodology. 
Hence we show next the results of applying the identification algorithm to actual measured 
data. However, even in the ideal low noise simulated cases, the true order of both closed loop 
systems could not be determined exactly although the frequency response of the identified 
systems matched quite accurately (< 1%). Of course this basic limitation is of concern to any 
identification method because internal variables or states generally do not contribute equally 
to the input /output data. The results demonstrate that this generic problem becomes more 
acute with higher levels of noise. In addition, the effects of nonlinearity can further blur the 
gap between dominant system states and those which appears as noise effects. 
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Figure 45: Input sensitivity, 100 minutes simulated data: True (dash), Identified 7th order (solid), 
high noise a = .03. 
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Figure 46: Singular values to determine order of closed loop transfer function across plant from 
simulated data, low noise a = 1CT 8 (left), high noise a = .03 (right). 
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Figure 47: Closed loop transfer function across plant, 20 minutes simulated data: True (dash), 
Identified 6th order (solid), low noise a = 10 -8 . 
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Figure 48: Closed loop transfer function across plant, 100 minutes simulated data: True (dash), 
Identified 7th order (solid), high noise a = .03. 
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A. 5 Identified models from experimental data 
A. 5.1 Identified model of (/ — KP)^ 1 

Figure 49 shows the time histories of commands, total inputs and the resulting outputs from 
a closed loop system identification experiment. Both command signals are designed to be 
wideband signals which are heavily weighted for frequencies below 5 Hz. This measurement 
data is used to compute the singular values which are used to determine a suitable order 
for the input sensitivity state space model. An upper bound order of 18 (same value used 
in simulation) is chosen to construct the data matrix which give these singular values. To 
increase the number of refinements by increasing the number of samples in the calculation of 
these singular values, an overlap of 2048 time points are used to produce 10 samples instead 
of 5 resulting in slightly improved results, as shown in Figure 50. The singular values show 
no significant gap to indicate a natural order of the system and a 9th order realization is 
chosen, guided by the 9th order of the analytical input sensitivity model given earlier. 

Figure 51 shows the frequency response of an identified 9th order input sensitivity model. 
The frequency response of this model roughly resembles the results from an analytical input 
sensitivity model over the frequency range of interest from .1 to 10 rad/sec. The figure 
show a typical dilemma in trying to decide with some confidence, which among the two 
approximate models is more accurate. 


A. 5. 2 Identifed model of (I — PK ) l P 

Figure 52 shows the singular values to determine an order of a state space model of the 
closed loop transfer function across plant. An upper bound order of 18 is chosen to compute 
these singular values and an 8th order realization is chosen. Again, there are no obvious 
gaps in the singular value spread but these singular values drops off more rapidly than the 
input sensitivity case. Hence we expect a slightly better fit of the data with 8 states. It is 
found that a 9tli order realization produced an unstable system and higher order realizations 
introduced unrealistic poles and zeros at unpredictable locations. 

Figure 53 shows the predicted and an identified 8th order model of the closed loop transfer 
function across plant. Note that this identified model only roughly matches the analytical 
frequency response. The match appears to be better in the second input, namely, paddle 
angle 5 P , than for the fan motor voltage input V m . This likely indicates a better analytical 
model in the dynamics associated with the paddle angle input, particularly in the prediction 
of the altitude z and angle of attack 6. 


A. 5. 3 Prediction error and summary 

We adopt the conventional view of an accurate model - one which accurately matches the 
measured output due to an arbitrary input, i.e. , “the essence of a model is its prediction 
aspect” (page 170, [5]). To this end, Figures 54 show comparisons of magnitude spectrums of 
measured, analytically predicted, and corresponding prediction errors. The prediction errors 
for analytical and identified 5 models shown in the figures are Z 2 norms over the frequency 
range .1 to 10 rad/sec normalized by the l 2 norms of the corresponding measured signals. 
Observe that the prediction errors for the identified model is less than that of the analytical 

5 The predicted response based on identified model used in the comparisons are actually two separate identified 
models obtained independently, input sensitivity and closed loop system across plant. However, this implies that 
a corresponding identified open loop plant which satisfies the closed loop equations cannot be determined exactly 
because it is overdetermined. 
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model but there is still significant discrepencies relative to the measured spectra. Note that 
the effects of significant but unknown disturbances and noise have not been included in the 
prediction by both models in the above comparisons. 

In summary, the identified closed loop model appears to be more accurate than the 
analytical model because it better predicts the measured signals. Supported by simulation 
results, it was found experimentally that the orders for the identified models for both input 
sensitivity and closed loop across plant were unclear because of the obvious presence of 
significant noise levels. Overall, the identified models resembles the analytical model but 
appears too large when this difference is viewed as an additive model error, for robust 
control application. 
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Figure 49: Closed loop measurements: command and total inputs (top), outputs (bottom) 
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Figure 51: Input sensitivity, 20 mins closed loop experimental data: Identified 9th order (solid), 
analytical (dash). 
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Figure 52: Singular values to determine order of closed loop transfer function across plant using 
experimental data, 10 samples each 4096 points at 20 Hz with 2048 overlap points. 
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Figure 53: Closed loop transfer function across plant, 20 minutes experimental data: Identified 8th 
order (solid), analytical (dash). 
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Figure 54: Comparison of measured and analytically predicted outputs and inputs 
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A. 6 Output noise model 


We consider a deterministic approach in the construction of an output noise filter from 
measured closed loop signals. The methodology is validated using simulated data before 
applying it to experimental data. Both low and high levels of simulated output noise are 
considered but the same command inputs are used in the simulation as in the experiments. 

Figure 23 illustrates the noise identification problem for a closed loop system. The inner 
loop controller, K is given (actually the same PD controller is used) and signals, y, u, 
and commands r are measured. The signal q represents an unknown equivalent additive 
output noise, “equivalent” because all unknown exogenous signals are modeled as noise at 
the plant output. This is a modeling assumption which is obviously not physically accurate 
but convenient in obtaining a mathematical model of the unknown exogenous signals in the 
system. Figure 23b shows the output noise effect as viewed from an open loop perspective. 

The goal is to obtain a suitable filter, V, such that q will be the output of V driven by a 
normalized discrete-time white noise, u. For simplicity, each channel of q is assumed to be 
independent so that V will be diagonal. To relax this assumption, further modeling effort 
is needed (for example the stochastic realization approach as described in [23] and [24]). 
As illustrated in Figure 23c, this noise model is the preferred form for subsequent model 
validation and structured uncertainty norm determination and subsequently to validate the 
performance of an outer loop controller, C. Notice that unlike the open loop case K = 0, 
output sensitivity is needed to properly determine the spectrum of q since simply fitting each 
closed loop output based on r = 0 could lead to significant errors when the output sensitivity 
matrix is not diagonally dominant. 

In this study, since we have no input test signal at the plant output to directly estimate 
the output sensitivity (so that one can effectively invert it to estimate q), we consider an 
approach which directly estimates the spectrum of q from frequency responses of the input 
sensitivity and closed loop across plant, which can be directly identified from available signals. 
A spectrum of q is obtained by solving the least squares problem corresponding to the 
overdetermined set of loop equations for inputs and outputs: 


where T ur := (/ — A'P) -1 , and T yr := (/ — PK)~ 1 P, are previously identified closed loop 
transfer functions obtained with independent experiments. Notice that solving for q uniquely 
from the output equation using an r = 0 experiment by inverting the output sensitivity 
matrix is actually a special case of the above least squares solution since / + T yr K, the 
coefficient matrix of q, is actually the output sensitivity matrix. Finally, to obtain V, an 
ensamble average of the absolute values of the least squares solution for each signal sample 
are computed and then fitted with a stable discrete filter for each channel. 

Due to feedback, the resulting closed loop output response due to the command input, r, 
and output noise, v, could be significantly different depending on the particular controller 
chosen, as evident from the equations 


In particular, the feedback will affect the signal-to-noise ratios which is a key parameter 
affecting the degree of success in system identification. We define the following signal-to- 
noise ratios at the outputs and inputs: 
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Notice that for zero mean signals, the above signal-to-noise ratios are simply the ratio of 
standard deviations of signals to noise. In principle, given input test signal, r, one can 
improve signal-to-noise ratios and hence the outcome of the closed loop system identification 
experiment by judiciously choosing the controller, K, such that the test signal response 
denoted by the first column partitioned transfer function matrix in equation A. 6 is maximized 
and the noise response denoted by the second column partitioned transfer function matrix 
is minimized, subject to closed loop stability during the test. 

Consider the following bounds on the signal-to-ratios: 
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The factor a (a) corresponds to the ratio of the minimum (maximum) output signal gain 
over maximum (minimum) noise gain at output, i.e. , the worst (best) signal to noise ratio at 
the output. Similarly, the factor (3 (J3) corresponds to the ratio of the minimum (maximum) 
input signal gain over maximum (minimum) noise gain at input, i.e., the worst (best) signal 
to noise ratio at the input. Figure 55 shows the signal-to-noise ratio lower and upper bounds 
a and a (solid) and /3, f3 (dash), which are based on the analytical plant model and a known 
controller. The bounds themselves are only a prediction to the extent that the plant model 
is an estimate of the physical system and that all unknown exogenous noise is represented by 
an equivalent additive noise model at the plant output. From the figure, notice that at low 
frequencies there is a wider range of possible variations in the signal-to-noise ratios and more 
so for the output channels. More importantly, note from the upper bounds of the signal-to- 
noise ratio predictions that at higher frequencies (say > 3 rad/sec), the best signal-to-noise 
ratio is less than 1, a bad sign for any sort of system identification or validation work based 
on measured signals at these frequencies. 


A. 6.1 Noise model from simulated data 

To validate the methodology used to realize a discrete filter driven by discrete-time white 
noise, simulated low and high output noise (similar to simulated noise used in system iden- 
tification study) was added internal to the feedback system consisting of PD controller and 
the analytical plant model. We assume that closed loop system identification has been done 
separately (see earlier simulation results in system identification) so that the identified in- 
put sensitivity and closed loop across plant are used, which do not match the truth model 
exactly. For the low noise case, the identified models (from simulated data) closely resemble 
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Figure 55: Lower and upper bounds of s/n ratios: a, a on output (solid) and f3, /3 on input (dash). 


the analytical closed loop system whereas in the high noise case, the identified closed loop 
models (from simulated data) differs substantially from the analytical model. 

Figures 56 show comparisons of identified output noise for low and high noise cases respec- 
tively. In all DFT computations, rectangular windows were used and the command input 
was zero in the closed loop simulation. At both noise levels, the noise magnitude spectrum 
was accurately recovered by this procedure at nearly all frequencies except at the very low 
frequencies where the conditioning of the least squares problem increases to 10 3 . Apparently 
this is a problem because the identified models were not very accurate at low frequencies 
which is amplified by the large condition number. 

In summary, the algorithm used for estimating an output spectra appears to work reason- 
ably well in simulation and so we apply it to actual measurement data. Using free response 
data (r = 0) also improved the recovery of the noise in the closed loop simulation. Appar- 
ently this avoids the Fourier Transform approximation error in the estimation of the effect 
of command input at the input and output by multiplying the identified frequency response 
matrices T yr and T ur by the DFT of r due to its time-limited signal. 


A. 6. 2 Noise model from experimental data 

An independent closed loop experiment with zero command input is carried out to capture 
the spectra of unknown exogenous noise in the closed loop system with the PD controller. 
A discrete time filter is fitted over the noise spectra for the purpose of model validation, 
uncertainty modeling, and designing a second controller over the PD control loop. The 
previously identified closed loop models are used in the least squares procedure for obtaining 
noise spectra as outlined earlier. Figure 24 shows the identified output noise spectra q (solid) 
and the corresponding stable rational filter fit for each output V (dash). To estimate the 
sensitivity of the least squares solution in obtaining the noise spectra, the condition number 
is also shown by the dash-dot line. Note the following observations: 
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• The condition number shows that the calculation in the noise spectra is less reliable 
at the lower frequencies for all three outputs. This sensitivity is with respect to both 
errors in the model and numerical. 

• Since the identified closed loop models used almost certainly is inaccurate, the noise 
estimate obtained from the least squares solution of the loop equations will include 
both noise and model error, i.e., the estimated noise is expected to be a conservative 
estimate. This implies a dependence of the least squares noise on the additive model 
uncertainty, and is consistent with the resemblcnce of the noise spectra to the frequency 
response of the identified model. 

• The noise spectrum show large disturbances at about 2.5 n rad/sec where n — 1, . . . and 
most prominent in the first output, (hjj/dtt and less in the third output, 6, and almost 
none in the second channel, z. This frequency matches the rotation rate of the Ducted 
fan about the z-axis. 

The following signal-to-noise ratios at the output and inputs are obtained: 

{s/n) y = (1.429,1.434,2.596) 

(s/n) u = (3.773,4.563) 

certainly, very noisy signals. 
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A. 7 Uncertainty models from simulated data 


Due to its inherent instability (see earlier description of Ducted fan in Appendix), any 
substantial experiment on the Ducted fan must be tested in closed loop. For the purposes of 
uncertainty bound determination and performance validation, for simplicity we postulate an 
uncertainty structure around the closed loop system with a PD controller as described earlier. 
Thus, the nominal model consists of the transfer function matrix from command input r to 
the closed loop output of the plant y, i.e. , (/ — PK)~ 1 P. The assumed independent additive 
equivalent output noise is amplified by the output sensitivity and filtered/normalized by the 
noise filter V. 


Again for simplicity, we choose unmodeled dynamics in both diagonal and unstructured 
forms of multiplicative uncertainties for all closed loop output channels. For this special case 
with no parametric uncertainties, zero disturbance allowance at input, and with only output 
multiplicative uncertainties, we obtain the following augmented plant: 



where G{P,K, V n ) 


0 [0, 0] (/-Pif)- 1 ? 

/ [0,U n ] ( I-PK)~ l P 


Since [Gn,Gi 2 ] = 0, the 77* (0,0 ) simplifies to 7/, = Gi 3 r, a constant for given data. The 
problem of finding a smallest model validating unmodeled dynamics then reduces to: 


min x 2 

subject to ||4(0, 4)|| 2 - x 2 |wi| 2 ||Gi 3 r|| 2 < 0, i = l,...,r 

m 2 < bi 


and since 4(0; 0) is affine and the inequalities represent ordinary norm bounds, a feasible set 
if it exists will be convex. For computational efficiency, we rewrite these as an optimization 
problem involving a linear cost function subject to a set of Linear Matrix Inequality (LMI) 
constraints [11]: 


min c T z 

Z 



subject to 

Qg ~ ||<4o,i|| 2 

Stz 

sym 

I 


\ b l 

sym 


Mz I 


where z := [Re(0); lm(0); Re(0); lm(0); x 2 ] G 'Jl 2n i>+ 2n <t>+ 1 and 


Qi ■= [ 

oj] 



Si ■ = 

Re(flj) — Im(flj) 

' Re (A) ' 

,0 

Im(Qj) Re(Oj) 

Im(A) 

M : = 

O2n0X2n^, I 2 ^212^x1 ] 


A := 

o 0 I n<j> jl n ^ 
j i n,p o o 



c ■ [O(2n^+2n0)i 1] 




a 2 := H 2 ||G 13 r|| 2 

The interior-point algorithm found in the MATLAB’s LMI Control Toolbox [20] is used to 
solve for the global minimum when feasible. 
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A. 7.1 Simulated data and case studies 


Figure 57 shows the true closed loop system and output noise which is used to generate sim- 
ulation data. This system consists of the unstable analytical plant model, Ptrue stabilized 
with a PD controller, K P D . The true noise, q. is simulated by a a sequence of zero mean 
Gaussian random number with standard deviation of a = 10” 8 for low noise and a = .03 for 
high noise. The command or test input r used in the simulation is chosen to be the same 
test input used in the experiments. 



r 


Figure 57: Truth model used for generating simulation data. 


Table 5 shows the different cases considered for simulating uncertainty bound determi- 
nation from data. To keep this study to a reasonable length, we consider only two forms of 
assumed uncertainty structure, namely, diagonal output multiplicative (denoted by “D” in 
Case labels) and unstructured or full output multiplicative (denoted by “F” in Case labels). 
These cases address, albeit in a limited sense, the effects of errors in the nominal model, 
uncertainty structure, noise allowances, and actual noise level in the data used for model 
validation. To enable us to define the correct or an incorrect uncertainty structure in this 
simulation study, the errors in the nominal model are simulated by perturbing an analytical 
truth model with speefied diagonal output multiplicative filters or with a full block output 
multiplicative filters, and by using identified models obtained through a system identification 
algorithm based on (low or high) noise corrupted signals. 

Cases Siml-1, Sim2-1, Sim3-1, Sim4-1 are meant to simulate a most ideal situation when 
the nominal model is a very accurate representation of the system of interest. Consequently, 
a very small model validating unmodeled dynamics is expected (if not zero) and serves as a 
necessary check of the model validation methodology. 

Cases Siml-2, Sim2-2, Sim3-2, Sim4-2 are meant to simulate the situation when the nom- 
inal model is inaccurate but has a correct uncertainty structure 6 . This is again an ideal 
situation since a correct uncertainty structure in a given application may not exist. Never- 
theless, these cases serve as a second necessary check of the model validation methodology 
wherein the correct uncertainty level is expected to be recovered. 

Cases Siml-3, Sim2-3, Sim3-3, Sim4-3 are meant to simulate the situation where the cor- 
rect uncertainty structure is unstructured. Clearly all cases that assume diagonal uncertainty 
structure (Sim*-3-D) will never recover the correct model error. Still, it does not rule out 
satisfying model validation conditions which involves only reproducing a particular output 
frequency spectra with the aid of an output noise allowance. 

In Cases Siml-4, Sim2-4, Sim3-4, Sim4-4, an identified model obtained from noisy identi- 
fication data is used as a nominal model and therefore there is no correct model uncertainty 

6 A philosophical argument can be made to the effect that in the modeling of physical systems, a particular 
uncertainty structure could be viewed as a particular parameterization of the discrepency between a mathematical 
model one uses and an imagined “true” model which one tries to contain using an uncertainty model. In other words, 
a “correct” uncertainty structure may not have any deep physical significance. 
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Case 

MV Data 

Nominal Model 

Uncertainty structure 

Noise allowance 

Siml-l-D 

Sim, low noise 

True 

Diagonal 

True peak 2-norm 

Siml-2-D 

Sim, low noise 

True + Diagonal error 

Diagonal 

True peak 2-norm 

Siml-3-D 

Sim, low noise 

True + Full block error 

Diagonal 

True peak 2-norm 

Siml-4-D 

Sim, low noise 

IDed (low noise data) 

Diagonal 

True peak 2-norm 

Sim2-1-D 

Sim, low noise 

True 

Diagonal 

True mean 

Sim2-2-D 

Sim, low noise 

True + Diagonal error 

Diagonal 

True mean 

Sim2-3-D 

Sim, low noise 

True + Full block error 

Diagonal 

True mean 

Sim2-4-D 

Sim, low noise 

IDed (low noise data) 

Diagonal 

True mean 

Sim3-1-D 

Sim, high noise 

True 

Diagonal 

True peak 2-norm 

Sim3-2-D 

Sim, high noise 

True + Diagonal error 

Diagonal 

True peak 2-norm 

Sim3-3-D 

Sim, high noise 

True + Full block error 

Diagonal 

True peak 2-norm 

Sim3-4-D 

Sim, high noise 

IDed (high noise data) 

Diagonal 

True peak 2-norm 

Sim4-1-D 

Sim, high noise 

True 

Diagonal 

.lx True mean 

Sim4-2-D 

Sim, high noise 

True + Diagonal error 

Diagonal 

.lx True mean 

Sim4-3-D 

Sim, high noise 

True + Full block error 

Diagonal 

.lx True mean 

Sim4-4-D 

Sim, high noise 

IDed (high noise data) 

Diagonal 

.lx True mean 

Siml-l-F 

Sim, low noise 

True 

Unstructured 

True peak 2-norm 

Siml-2-F 

Sim, low noise 

True + Diagonal error 

Unstructured 

True peak 2-norm 

Siml-3-F 

Sim, low noise 

True + Full block error 

Unstructured 

True peak 2-norm 

Siml-4-F 

Sim, low noise 

IDed (low noise data) 

Unstructured 

True peak 2-norm 

Sim2-1-F 

Sim, low noise 

True 

Unstructured 

True mean 

Sim2-2-F 

Sim, low noise 

True + Diagonal error 

Unstructured 

True mean 

Sim2-3-F 

Sim, low noise 

True + Full block error 

Unstructured 

True mean 

Sim2-4-F 

Sim, low noise 

IDed (low noise data) 

Unstructured 

True mean 

Sim3-1-F 

Sim, high noise 

True 

Unstructured 

True peak 2-norm 

Sim3-2-F 

Sim, high noise 

True + Diagonal error 

Unstructured 

True peak 2-norm 

Sim3-3-F 

Sim, high noise 

True + Full block error 

Unstructured 

True peak 2-norm 

Sim3-4-F 

Sim, high noise 

IDed (high noise data) 

Unstructured 

True peak 2-norm 

Sim4-1-F 

Sim, high noise 

True 

Unstructured 

.lx True mean 

Sim4-2-F 

Sim, high noise 

True + Diagonal error 

Unstructured 

.lx True mean 

Sim4-3-F 

Sim, high noise 

True + Full block error 

Unstructured 

.lx True mean 

Sim4-4-F 

Sim, high noise 

IDed (high noise data) 

Unstructured 

.lx True mean 


Table 5: Simulated cases: effects of nominal model, uncertainty structure, and noise allowance. 
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structure. This ambiguity is whats expected in a typical application. This contrasts sharply 
to the other cases where a correct uncertainty structure exists, albeit simulated, and the 
assumed uncertainty structure is either correct or incorrect. Still, analogous to the previous 
case, the existence of a model validating set cannot be ruled out, even if a “correct” model 
uncertainty structure is unknown. However, the identified model from simulated data closely 
matches the true model. 


A. 7. 2 Perfect nominal model 

Figure 58 shows diagonal uncertainty bounds for cases Siml-l-D,. . ., Sim4-1-D which assumes 
perfect nominal model. This is a reference case where the recovered uncertainty bounds 
should be very small or ideally zero. With a low level of noise in the data, cases Siml-l-D 
and Sim2-1-D show a very small level of unmodeled dynamics (circle and square dash lines). 
This very small level of unmodeled dynamics obtained is independent of the noise allowance 
levels assumed, since the circle and square dash lines overlap. However, notice that even 
with a perfect nominal model and the data contains only a small level of noise, the model 
validating unmodeled dynamics is not exactly zero. This residual is due to the fact that the 
output error used in model validation definition assumes 

e y {z) := Z[y(t)\ - T yr (z)Z[r(t)\ = 0 

under the ideal conditions: perfect model T yr (z), no noise or disturbance effects in y{t) and 
r, and signals are either infinite time length or is periodic. However, since all measured time 
signals are of finite length even under near ideal conditions where noise effects are negligible, 
the following output error results 

e y {z) := Z[y(t)h(t )] - T yr (z)Z[r{t)h(t)\ = Z[y(t )] * Z[h(t )] - T yr (z ) (. Z[r(t )] * Z[h(t)}) % 0 

where Z\-\ denotes the z-transform operation, * denotes the frequency convolution operation 
and h(t) denotes the chosen time window (refer earlier section for a more detail discussion on 
the assumptions on signals and systems). Figure 59 illustrates this residual error for a perfect 
model with a Hanning window over a 12000 point (600 seconds) sample, which is exactly 
half the total length of experimental data. In this application, this residual DFT error is 
considered an acceptably small level of fictitious model validating unmodeled dynamics. 

Figure 58 also shows what happens to the uncertainty bounds when there is a high level 
of noise in the data used for model validation. Using a high noise allowance level equal to the 
peak true, Case Sim3-1-D shows a very small uncertainty bound (diamond dot line), which 
is acceptable. Apparently, the remaining noise allowance is sufficient to compensate for the 
DFT errors. However, if only a small level of noise allowance (10% of true mean noise) is 
used when there is a high level of noise in the data, as in case Sim4-1-D, a large erronous 
level of uncertainty bound results (triangle dot line), mostly to account for the large output 
noise. 

Assuming an unstructured uncertainty with a perfect nominal model, the computed un- 
structured uncertainty bounds obtained are similarly very small as shown by the circle and 
square dotted lines in Figure 60. The erroneous effects of using a low level of noise allowance 
when there is a large level of noise in the data is similar (triangle dot line). 


A. 7. 3 Nominal model with correct uncertainty structure 

Figure 61 shows the diagonal uncertainty bounds for cases Siml-2-D, Sim2-2-D, Sim3-2- 
D, Sim4-2-D, whose nominal models are chosen as the analytical truth model perturbed 
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diagonally at the output as follows 


(37) 


P{z)nom = {I + A (z))P true (z) 

where 


A(z) : = 

*) ■■= 

A diag ■ 

This means that a correct uncertainty structure and level exists. The solid lines in Figure 
61 shows the correct model error levels, a constant in the first channel and first order low 
and high pass filters for the second and third channels, respectively. The results show that 
with low noise levels in the data as in Siml-2-D and Sim2-2-D, the correct diagonal model 
error levels are recovered reasonably accurately as indicated by circle and square dash lines 
in Figure 61. At most frequencies, the predicted uncertainty bound is slightly less than the 
correct levels (except at the very low frequencies where the data record length becomes a 
limiting factor) which do not appear to be due to the noise allowance levels (since the circle 
and square dash lines overlap). Figure 61 also shows that if a high level of noise is present 
in the data for model validation, as in cases Sim3-2-D and Sim4-2-D, apriori knowledge of 
the correct uncertainty structure does not help to recover the correct uncertainty bounds 
(compare line with diamond and triange dash lines). 

As a second example of a situation where a correct uncertainty structure is known, Fig- 
ure 62 shows the computed unstructured uncertainty bounds for a nominal model with a 
full block output model error and assuming correctly, an unstructured uncertainty at the 
output. For the cases with low noise in the data, Siml-3-F and Sim2-3-F, the uncertainty 
bounds (circle and square lines) roughly recovers the maximum singular value of the correct 
unstructured model error (solid line). However, with large levels of noise in the data as in 
cases Sim3-3-F and Sim4-3-F, erroneous levels of uncertainty bounds result (diamond and 
triangle lines) independent of the level of noise allowance assumed. 



W{z) A diag 

dia 4 2 - H-S 

-^3x3) ^"(Arfjag) 1 


.5208 


(z - .7776) 
(z + .1584) 


A. 7. 4 Nominal model with incorrect uncertainty structure 

Figure 63 shows the diagonal uncertainty bounds for cases Siml-3-D, Sim2-3-D, Sim3-3-D, 
and Sim4-3-D, whose nominal model is the truth model perturbed at the output as follows 

P{z) nom = (/ + A(z))P true (z) 


where 






A(z) 

:= W(z) A full 






' -0.0896 

- 0.2473* 

—0.0178 + 0.5314* 

0.4692 - 0.2975* 

< 

: = 

-0.1352 

- 0.1399* 

-0.3587 + 0.1471* 

-0.1014- 0.2886* 


0.0674 + 0.4169* 

-0.1886 - 0.2629* 

-0.1121 + 0.2525* 


o’(Afuii) — 1 


and W(z) is a diagonal filter denoting the correct uncertainty radius for each output. This 
is an example of a situation where the assumed uncertainty structure is incorrect and is 
incapable of recovering the correct uncertainty bounds. Figure 63 shows the uncertainty 
bounds obtained for all four cases which do not resemble the maximum singular value of the 
correct unstructured output model error indicated by the solid line. Not surprisingly, due to 
the nonuniqueness of model validating uncertainties, all four cases satisfied model validating 
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conditions in spite of the incorrect uncertainty structure assumed. Finally, note that for 
the two cases with low level of noise in the data (Siml-3-D and Sim2-3-D), the uncertainty 
bounds are very insensitive to noise allowance (circle lines overlap square lines). 

As an example of a situation where the assumed uncertainty structure is incorrect because 
it is overparameterized, Figure 64 show the results of Siml-2-F,..., Sim4-2-F, where the 
nominal model is diagonally perturbed as given by equation 37, which is identical to cases 
Siml-2-D,..., Sim4-2-D. For the case with low levels of noise in the data (Siml-2-F and Sim2- 
2-F), the maximum singular value of the unstructured uncertainty (circle and square dash 
lines in Figure 64) roughly matches the maximum singular value of the correct diagonal 
model error (solid line). Of course this successful recovery is likely due to the fact that 
the assumed unstructured (and overparameterized) uncertainty is capable of generating the 
correct diagonal uncertainty bounds. Again, with high levels of noise in the data (Sim3-2-F 
and Sim3-2-F), erroneous uncertainty bounds result (diamond and triangle lines in Figure 
64). 


A. 7. 5 Identified nominal models 

In this set of cases, identified models are used as nominal models. Cases (Siml-4-D, Sim2-4- 
D) and (Siml-4-F, Sim2-4-F) uses nominal models obtained from system identification based 
on signals which contain low levels of noise, hence these are more accurate identified models. 
The complementary cases (Sim3-4-D, Sim4-4-D) and (Sim3-4-F, Sim4-4-F) uses inaccurate 
identified nominal model due to the use of system identification data containing a high level 
of noise (compare model errors in Figures 47 and 48). 

Given a 3 by 2 truth model, Pt rU ei and an identified nominal model P id , a correct additive 
uncertainty can be defined and computed unambiguously, but a multiplicative output uncer- 
tainty can only be defined by solving for A mu i t from the following underdetermined matrix 
equation: 

(-^3x3 + A )P id Ptrue 

Notice that there is no physical preference for a least squares A corresponding to an output 
multiplicative uncertainty (let alone correct uncertainty structure). This is an example where 
we are only certain that the two transfer matrices are different and the correct uncertainty 
structure is therefore indeterminate. 

Figure 65 shows the unstructured output multiplicative uncertainty bounds (circle and 
square dash lines) for cases Siml-4-F and Sim2-4-F which uses a low level of noise in data. 
Figure 66 shows the diagonal output multiplicative uncertainty bounds (circle and square 
dash lines) for cases Siml-4-D and Sim2-4-D which uses a low level of noise in data. Both 
diagonal and unstructured output multiplicative uncertainty bounds compares surprisingly 
well to the maximum singular value of the least squares output multiplicative uncertainty 
as shown by solid lines in Figures 65 and 66. In spite of the unknown correct uncertainty 
structure for the identified models, all four cases satisfied model validation conditions and 
produced reasonable uncertainty levels. 

Figures 67 and 68 show the uncertainty bounds using inaccurate identified nominal mod- 
els and data with high levels of noise. Clearly the uncertainty bounds (diamond and triangle 
dotted lines) do not resemble the maximum singular value of the least squares output mul- 
tiplicative uncertainty as shown by solid lines. Again, the predicted uncertainty bounds are 
sensitive to the assumed noise allowance because of the high levels of noise in data. 
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A. 7. 6 Summary 

Since an infinity of other uncertainty structures, parameters, noise allowance levels, and 
actual noise in data are possible but cannot all be considered in any simulation, the results 
are in no way conclusive. Nevertheless, the results suggest the following: 

For all following cases with a low level of noise in data, all uncertainty bounds were 
insensitive to the noise allowances used. With accurate nominal models, the uncertainty 
structure was not important because the resulting uncertainty bounds required were as ex- 
pected small. Namely, for the perfect nominal case, assuming diagonal uncertainty structure 
resulted in very small levels (circle and square lines in Figure 58) while an unstructured un- 
certainty structure assumption also resulted in very small levels of maximum singular value 
uncertainty as expected (circle and square lines in Figure 60). For the accurately identified 
model case, assuming unstructured uncertainty gave matching small levels of maximum sin- 
gular value uncertainty levels (circle and square lines in Figure 65) while assuming diagonal 
uncertainty structure resulted in matching small levels of diagonal uncertainty levels (circle 
and square lines in Figure 66) compared to maximum singular value of the least squares un- 
structured output multiplicative uncertainty (solid line in Figure 66). These fictitious small 
levels of uncertainty bounds consists mostly of DFT errors and were mitigated by judiciously 
windowing reasonably long samples. 

For less accurate nominal models with correct a priori assumptions on the uncertainty 
structure, the correct uncertainty levels are recovered for both diagonal (circle and square 
lines compared to solid line in Figure 61) and unstructured cases (circle and square lines 
compared to solid line in Figure 62). The underestimated uncertainty levels at most fre- 
quencies are due to the fact that the particular sample data did not reflect the worst case 
directional response for this multivariable system. 

For less accurate nominal models with incorrect a priori assumptions on the uncertainty 
structure, the uncertainty bounds obtained were generally unpredictable. This unpredictabil- 
ity is of course not totally surprising since the incorrect uncertainty structure will not allow 
the recovery of the correct uncertainty bounds. For example, in the nominal model case with 
unstructured error at output, assuming diagonal uncertainty led to unrecognizable bounds 
(circle and square lines compared to solid line in Figure 63). However, notice from the fig- 
ure that at each frequency, the maximum uncertainty component (maximum singular value 
of a diagonal matrix) apparently is greater than the maximum singular value of the cor- 
rect full uncertainty matrix. As a second example with incorrect a priori assumptions on 
the uncertainty structure, consider the overparameterized uncertainty case where the correct 
uncertainty is diagonal but the assumed uncertainty model is assumed unstructured. The 
maximum singular value of the uncertainty bounds of the assumed unstructured uncertainty 
(circle and square lines in Figure 64) compares reasonably well, albeit an underestimation, to 
the correct maximum uncertainty component (maximum singular value of the true diagonal 
marix) indicated by solid line. This somewhat successful recovery of the correct uncertainty 
bound ' was made possible because the true diagonal uncertainty structure is a subset of 
the corresponding unstructured uncertainty. Interestingly, note from the two previous ex- 
amples that the maximum simgular value of the unstructured uncertainty was consistently 
smaller than that of the diagonal uncertainties, which suggests that assuming less structure 
in the uncertainty, and therefore intuitively more conservative, may not necessarily lead to 
a more conservative uncertainty model, by reducing uncertainty values with more uncertain 
variables. 

For black-box type nominal models typical of identified models, a correct uncertainty 
structure is indeterminate, given only the identified noiminal model and possibly an imag- 
ined true model. Due to this indeterminacy, any uncertainy bound comparisons are suspect. 

' although apples and oranges are compared since a diagonal matrix is compared to a full matrix 
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Nevertheless for this study, we computed and compared an equivalent output multiplicative 
error based on a minimum norm least squares solution to the output multiplicative uncer- 
tainty bound from data. Surprisingly, both the unstructured uncertainty bound in (circle 
and square lines in Figure 65) and diagonal structure uncertainty bounds (circle and square 
lines in Figure 66), match reasonably well with the maximum singular values of the least 
squares equivalent output multiplicative uncertainties (solid lines in Figure 65 and 66). In 
summary, incorrect or indeterminate uncertainty structure did not preclude satisfying model 
validation conditions but sometimes resulted in unpredictable uncertainty norm bounds. 

For all cases with a high level of noise in data , the uncertainty bounds were unreliable and 
sensitive to the assumed levels of noise allowance. In the selection of noise allowance levels, a 
larger level is preferred over smaller levels (than true noise level in data) to avoid erroneous 
uncertainty bounds, unless the noise in the data is known to be extremely low. Although 
it is intuitively clear that any attempt to validate a model will be less reliable with higher 
noise levels, as shown here, the results have been further subject to a limiting feature of the 
model validation theory investigated which do not distinguish noise and model error effects. 
In fact, the validation approach is to use noise variables as a given allowance to minimize 
the model error uncertainty bounds necessary to satisfy model validation conditons. Hence 
further research is recommended whereby the noise signals in the data can be discriminated 
from model error effects in satisfying model validation conditions so that noise allowance 
will not compensate for model errors and vice versa. 
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Siml —1 — F(cir — ),Sim2-1-F(sq — ),Sim3-1-F(dia:), Sim4— 1 — F(tri:), han, 12000 pts 



Frequency (rad/sec) 


Figure 60: Unstructured uncertainty bound based on perfect nominal model, Siml-l-F,...,Sim4-l-F. 
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Sim 1 — 3— F(cir — ), Sim2-3-F(sq — ), Sim3-3-F(dia:), Sim4— 3— F(tri:), MSV True error(line), han, 12000 pts 



Frequency (rad/sec) 

Figure 62: Unstructured uncertainty bound based on nominal model with full block error, Siml-3- 
F,...,Sim4-3-F. 
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Sim 1 — 2— F(cir — ), Sim2-2-F(sq — ), Sim3-2-F(dia:), Sim4— 2— F(tri:), MSV True error(line), han, 12000 pts 



Figure 64: Unstructured uncertainty bound based on nominal model with diagonal block error, 
Siml-2-F,...,Sim4-2-F. 


Sim1-4-F(cir — ), Sim2-4-F(sq — ), max SV of Least Sq (line), han, 12000 pts 



Figure 65: Unstructured uncertainty bounds for identified nominal model under low noise: Sirnl- 
4-F and Sim2-4-F. 
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Sim3— 4— F(dia:), Sim4— 4— F(tri:), max SV of Least Sq (line), han, 12000 pts 



Figure 67: Unstructured uncertainty bounds for identified nominal model under high noise: Sim3- 
4-F and Sim4-4-F. 
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